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Abstract 

We present an overview of the evolution of ab initio methods for few-nucleon systems with 
A> A, tracing the progress made that today allows precision calculations for these systems. First 
a succinct description of the diverse approaches is given. In order to identify analogies and dif- 
ferences the methods are grouped according to different formulations of the quantum mechanical 
many-body problem. Various significant applications from the past and present are described. 
We discuss the results with emphasis on the developments following the original implementations 
of the approaches. In particular we highlight benchmark results which represent important mile- 
stones towards setting an ever growing standard for theoretical calculations. This is relevant for 
meaningful comparisons with experimental data. Such comparisons may reveal whether a specific 
force model is appropriate for the description of nuclear dynamics. 



> : 1 Introduction 



5-h ■ The last two decades have witnessed decisive progress in nuclear physics by ab initio calculations. In 



the previous years ab initio calculations were restricted mainly to applications in two- and three-nucleon 
problems while in the recent past these techniques have been extended to considerably larger nuclei, 
which, in some cases, consist of even more than ten nucleons. Here we include such larger systems in 
our definition of few-nucleon systems. Progress has been made in various directions. Methods which 
were already well established in the field have been further improved. As examples we mention Faddeev 
calculations, which have been extended from three- to four-nucleon systems [1, 2] and which have solved 
the problem of inclusion of the Coulomb interaction in momentum space [3] , and Green's function Monte 
Carlo (GFMC) calculations, which have reached a degree of sophistication which allows inclusion of local 
realistic nuclear forces up to a nucleon number of A = 12 [4]. In addition, novel approaches such as 
the no-core shell model (NCSM) [5] and the Lorentz integral transform (LIT) [6] methods have been 
introduced. Last, but not least, growing computational resources and improved numerical algorithms 
have allowed calculations which were previously unthinkable. 

At the same time enormous progress has also been made in the description of the dynamical input of 
an ab initio calculation, i.e. the nuclear interaction. In the 1990s modern realistic phenomenological and 
meson-theoretical nucleon-nucleon (NN) potential models, which provide high precision descriptions of 
both NN bound state and the wealth of NN scattering state observables, have been developed. However, 
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the failure of these potentials to account for the binding energy and other observables in the three-body 
sector has demonstrated the necessity of including three- and in principle even more-nucleon forces. 
Three-body forces of different forms have been proposed and their parameters have been obtained by 
fitting to data. At the same time a new line of research dedicated to the derivation of the nuclear 
interaction from chiral perturbation theory was opened up. This new ansatz has led to the construction 
of competitive modern realistic potentials that contain chiral many-nucleon forces consistent with the 
two-body force [7, 8, 9]. 

The complications in ab initio calculations, generated on the one hand by the strongly repulsive 
nature of all the modern realistic potentials of the 1990s and on the other hand by the three-body forces, 
have stimulated an interesting discussion and led to a corresponding line of research. Thus the last few 
years saw considerable evolution in obtaining phase-shift equivalent NN potential models, which often 
lead to a great simplification in ab initio calculations. In fact one may apply unitary transformations 
to the two-body Hamiltonian, and to other operators at the same time, without effecting any change 
in the two-nucleon observables. This leads to an infinite number of equally valid NN potentials some 
of which could have characteristics that are simpler to treat in numerical calculations. For example 
one can aim at softer potentials at short range. By applying a unitary transformation at the three- 
body level one could even trade a three-nucleon force with a two-body potential that has a different 
off-shell behavior. Examples for this procedure are the potentials obtained by the renormalization 
group transformations like the so-called V\ ow k [10], and the free or in medium similarity renormalization 
group (SRG) potentials [11], as well as the force obtained by the unitary correlation operator method 
(UCOM) [12] and the J- matrix inverse scattering potential (JISP) [13]. 

Though it is true that all the various potentials mentioned above are phase-shift equivalent, they 
may not be identical with the original potential models, since, in general, the off-shell structure can be 
different. This may lead to different results for observables in the A > 2 sector, such as e.g. the triton 
binding energy. It is clear that even applying a specific unitary transformation at the three-body level, 
one may have different results for A > 3 observables and one may regard the transformed Hamiltonian 
as a new potential model. Since there are an infinite number of different unitary transformations one has 
an infinite number of different potentials. Only the comparison to experiment of precise ab initio results 
for A > 3 observables can discriminate among these potentials. In this work we will not discuss nuclear 
forces, although in various places we do mention different models for NN potentials and three-nucleon 
forces (for a recent review of the latter see [14]). 

The present work aims at giving an overview of the, by now, numerous ab initio methods which are 
applied to obtain ab initio solutions of nuclear many-body problems for bound and scattering states, 
as well as for observables of hadronic and perturbation-induced reactions. Since in the last years the 
term ab initio has been rather loosely used, we first would like to clarify what we consider to be an 
ab initio method and an ab initio result. For this purpose let us consider an A nucleon system that is 
described by a well-defined microscopic Hamiltonian H with A nucleon degrees of freedom and where 
the internal relative motion is treated correctly. If a method enables one to obtain the observable 
under consideration by solving the relevant quantum mechanical many-body equations, without any 
uncontrolled approximation, we consider it to be an ab initio method. Controlled approximations, 
however, are allowed. In fact a controlled approximation, e.g. a limited number of channels in a 
Faddeev calculation, can be increasingly improved up to the point that convergence is reached for the 
observable. Such a converged result we denote as a precise ab initio result. The comparison of precise 
ab initio results with nuclear data then allows an indisputable answer as to whether or not the chosen 
Hamiltonian appropriately describes the nuclear dynamics. Any uncontrolled approximation in the 
calculation would not lead to such a clear-cut conclusion. Quite naturally, precise ab initio results 
obtained with different ab initio methods but with the same Hamiltonian as input, have to agree and 
are often referred to as benchmark results. 

The aim of a quantum mechanical treatment of a particle system is the determination of observables. 



Often this aim can be achieved by rather different formulations. The organization of our work is guided 
by this diversity. In the present few-body case all formulations start from the same point, namely from 
a Hamiltonian that has nucleons as degrees of freedom and that contains a well-defined interaction 
model. While pursuing a specific formulation one may still have a choice between different methods, 
some of which can be innate to the strategy under consideration, or others just adapted from elsewhere. 
Therefore the first part of the work is devoted to such different formulations. The diversity of methods 
are described in the corresponding subsections, which contain, in addition, some illustrative examples 
from the literature. Selected results are presented and discussed in a following section. 

Here a comment is in order on the issue of whether a given method is more suited to bound- and/or 
continuum-state problems. In principle all the various methods presented here can be applied to both 
cases. The only exceptions are the LIT and the complex scaling methods, which do not treat bound- 
state problems, but which are dedicated to the calculation of observables in the continuum. On the 
other hand it is true that some methods are more often used, either for bound-state or for continuum- 
state problems. For instance, variational techniques are applied more frequently to bound states, but, as 
discussed in Section 2.3, continuum-state calculations can also be carried out via the Kohn variational 
principle. A similar observation can be made for the GFMC technique (see Section 2.5.1). 

The methods we are about to present in this work are quite general and not restricted to a specific 
field. However, we limit our discussion mainly to examples from nuclear physics. In order to make the 
wide range of applications more explicit we list here topics and related review articles, where ab initio 
techniques are receiving a growing interest from the few- and many-body communities: universality 
in few-body systems and Efimov physics [15], few-body ultracold atomic and molecular systems [16], 
quantum degenerate Fermi gases [17], and last but not least quantum dots [18]. For the latter case there 
is also a recent benchmark calculation with diffusion Monte Carlo and coupled cluster techniques [19]), 
where in addition a brief overview on the state-of-the-art of the physics of quantum dots is given in the 
introduction. 

In more detail, our work is organized as follows. In Section 2 we characterize the so-called few-body 
problem very briefly. In Section 2.1 we start with short introductions to the Faddeev and Faddeev- 
Yakubovsky formalisms. A reformulation of the Faddeev approach and the Alt-Grassberger-Sandhas 
(AGS) method are discussed in Section 2.2. An introduction to the variational formulation of the 
few-body problem is given in Section 2.3, followed by various subsections describing specific variational 
methods. Another formulation, largely used both in few-and many-body physics, especially for the 
construction of effective interactions, is based on similarity transformations. Thus, in Section 2.4, 
the similarity transformation formulation is outlined and few-body methods that make use of it are 
described. Formulations that are suitable for using Monte Carlo techniques are discussed in Section 2.5. 
After the discussion of the various ab initio methods, we show in Section 3 a selection of results that in 
our view represent some advancement of the field. There it will be seen that we place a certain emphasis 
on benchmark calculations. In Section 3.1 we report on the early stage of ab initio calculations for the 
A = 4 nuclear system, because there the fundamentals for many later developments have been laid and 
we think that this pioneering work deserves not to be forgotten. The more modern stage of A = 4 
calculations initiated in the late 1980s is discussed in Section 3.2. Results for nuclear systems with 
A > 4 are described in Sections 3.3 and 3.4. Finally summary and outlook are given in Section 4. 

2 The Few-Body Problem 

The dynamics of a system of A non-relativistic nucleons of equal mass m is governed by the translation 
invariant nuclear Hamiltonian H, which consists of kinetic energy T and potential V: 

H = T + V 
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A 2 t>2 A A 

In the equation above is the momentum of the i-th nucleon, Pcm is the center of mass (CM) 
momentum, while V^- and denote the NN potential Vnn, and the three-nucleon potential Vnnn, 
respectively. We will refer to the latter as three-nucleon force (3NF). Since the nucleons are effective 
degrees of freedom, the nuclear Hamiltonian H should in principle contain potentials involving more 
than three nucleons. However, at present such many-nucleon forces seem to play a negligible or at least 
a minor role in the dynamics of A- nucleon systems (see e.g. [7]). 

The basic dynamical equation to be solved ab initio is the Schrodinger equation, 

(H-E)\V) = 0, (2) 

where E and denote one of the eigenenergies and eigenfunctions of H, respectively. The spectrum 
of H, represented by the infinite set of eigenenergies E, is discrete below, and continuous above, the 
first break-up threshold (-Eth) of the A-nucleon system. 

In order to solve the Schrodinger equation one has to supply proper boundary conditions. For 
example, in coordinate space calculations one requires the regularity condition at = 0, where is 
the relative distance between two nucleons. The asymptotic boundary condition depends on the energy 
of the system. For E < E t ^ the wave function represents a bound state and thus is described by a square 
integrable (localized) function. For E > E t ^ the asymptotic boundary condition for the corresponding 
scattering state remains quite simple in the case of A = 2. On the contrary the asymptotic boundary 
conditions pose a serious problem if both E > E th and A > 2. How one solves that problem is discussed 
in the following subsection. In the A = 2 case one may start not from the Schrodinger equation but 
from its integral form, the Lippmann-Schwinger (LS) equation 

\*&) = \*) + E± l_ T V\*W), (3) 

where |$) is the free solution describing a plane wave, i.e. satisfying (T — E)\$) = 0. The nice feature 
of the LS equation is that it already contains the asymptotic boundary condition, i.e. an outgoing and 
an incoming spherical scattered wave for |^ (+) ) and |# (_) ), respectively [20]. 

This simple example shows that it can be helpful to reformulate the original quantum mechanical 
many-body problem so that it becomes possible to tackle it with a different technique. It is this 
possibility that generates the richness of methods in many-body theory. 



2.1 The Faddeev and Faddeev-Yakubovsky (FY) Formulation 

Today it is well known that a direct application of LS-type equations to a scattering problem does 
not lead to a unique solution for a system with more than two particles. However, how such a unique 
solution can be obtained was not clear for quite some time. In his seminal work Faddeev showed how 
the problem can be treated [21]. Considering the three-body system he derived a set of equations, the 
Faddeev equations. We do not illustrate Faddeev's original derivation, but discuss here a few selected 
points of the didactic derivation of the Faddeev equations given by Glockle [22]. The important result 
of Faddeev is that a specific three-body scattering state should obey not just one LS equation, but two 
additional equations at the same time. Rewriting the Hamiltonian as 

H^T + j^Vi, (4) 

i=i 
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where V 1 = V 23 , V 2 = Vi 3 , V 3 = Vi 2 , the three equations for the (2 + 1) — > (2 + 1) reaction that ensure 
the correct boundary conditions read as follows 



l*i +) > = |*«> + G a (E)V a \^) , (5) 
l*i +) > = Gp(E)V^M), (6) 

|*«> = G 7 (£0^|*W> , (7) 

where a 7^ /3 7^ 7 can take the values 1, 2, and 3 , is the total energy of the three-nucleon (3N) system 
and |^ + ^) is the outgoing scattering wave in channel a (i.e. where asymptotically one has particle a 
and the pair formed by the other two particles). This function satisfies 

#l*i +) > = £|*i +) >, (8) 
while is eigenfunction of the channel Hamiltonian H a = T + V a , i.e. it fulfills 

H a \$ a ) =E\$ a ). (9) 
The potential term V a is defined as V a = H — H a and G a (z) is the resolvent operator for channel a 

Ga(z) = ~ -TF" • (10) 

-2 — tl n 



The inhomogeneous LS equation (5) alone does not define the boundary conditions uniquely. This is the 
case only with the two additional homogeneous equations (6) and (7). Also for the (1 + 1+1) scattering 
one arrives at three analogous equations. However, all these sets of equations do not directly allow 
practical solutions. In fact for any equation set the same state has to fulfill three different equations. In 
order to come to a more feasible task the three separate equations should be transformed into a coupled 
equation system. The proper strategy is to decompose the wave function in three parts, the so-called 
Faddeev components \^ a ^): 

l*i +) > = EIVw> , (ii) 

where \4> a ,n) = Go(E)Vn\'if^) with G representing the free propagator. Using this decomposition one 
arrives at the following coupled equation system 

\* a ,a) = \<5> a ) + G a (E)V a (\y a ,p) + \V a ^)), (12) 
|* aJJ > = G p (E)ve(\* a , 7 ) + \* a , a )), (13) 

|tt a>7 > = G 7 (E)V^\^ a , a ) + |* a>/J » . (14) 

These are then the equations that are usually called Faddeev equations. For bound states one has a 
modified set of Faddeev equations, where the inhomogeneous term in Eq. (12) is dropped. 

The Faddeev equations in the form discussed above are naturally suitable for a solution in momentum 
space. In fact most applications for the three-nucleon (3N) system have been carried out in this 
representation (for a review see [23]). Also Faddeev calculations for the four-nucleon (4N) bound state 
have been made in momentum space as discussed below. Coordinate space Faddeev calculations have 
been performed for both the 3N and 4N cases. A particularly interesting example of using the Faddeev 
formalism in coordinate space is in low-energy (3+1) scattering. Therefore we devote the rest of this 
section to the Faddeev formalism in coordinate space. First we consider the 3N system and assume 
that the nucleons interact via a short-range two-body potential V^. 

Denoting with the position of nucleon i one has the following sets of Jacobi coordinates 

■Xij = Tj - Vi , y ij>k = r k - . (15) 



There are six different permutations and hence six different sets. However, x^- and x 3 -j describe the same 
physical situation, thus there are only three independent sets. In the following step the wave function 
is decomposed into the aforementioned three Faddeev components ipf 

(xy \^) = ^i(X23, Y23,l) + ^2(X 3 1, y 3 l,2) + ^3(Xi2, yi2,3) , (16) 

where we have suppressed spin and isospin degrees of freedom. Then, one requires that the Faddeev 
components satisfy the following set of coupled equations, 

(T-£ + \/ 2 3)^l(x23,y23,l) = -^23 (^2(x31,y 3 l, 2 ) +^3(Xl2,yi2,3)) , (17) 
(T- E + V 13 )ljj 2 (X3i,y 3 i i2 ) = -^13(V , l(x23,y23,l)+^3( x 12,yi2,3)), (18) 
(T-E + Vi 2 )^ 3 (Xl2,yi2,3) = -^12 (^l(x 23 ,y23,l) +^2(x 3 l,y 3 l,2)) , (19) 

which are also known as Faddeev-Noyes equations. Adding up the three Eqs. (17-19) one easily verifies 
that the Schrodinger equation, 

{T-E + V 23 + V 13 + V 12 ){x.y\V) = 0, (20) 

is fulfilled. For the case of identical particles (the nucleons) considered here, there is a nice feature 
of the Faddeev formalism, namely that the various ipi can be obtained from each other by particle 
permutations. In fact one has 

^2(x 3 l,y31,2) = ^12^23^1 (X 2 3,y 23, l) , ^z{*12, Yvifi) = ^23^12^1 (x 2 3, y23,l) , (21) 

and therefore it is sufficient to solve a single equation 

(T-E + V 23 )M*23, y 23 ,i) = ~V 23 (^12^23 + P23P12) ^i(x 23 , y 23 ,i) ■ (22) 

Let us consider a 3N state with total angular momentum J and parity n. To solve Eq. (22) for this 
state | J 71 ") the Faddeev component ipi, from here on simply denoted by tp, is expanded into channels 

^,y) = E^ M / w (^) ) (23) 

{A} X V 

where for simplicity the subscripts of variables x and y and the dependence of the functions ip and 
/1 A 1 on spin and isospin variables are dropped. Here the different channels refer to different sets 
{A} of quantum numbers. As to boundary conditions, one has (f)^(x =0, y) = (p^ x \x, y = 0)=0 for 
bound and scattering states. Moreover, in the bound-state case, one has the asymptotic boundary 
condition <f)^(x — X mscK , y)=(p^(x, y = Y mSLX ) =0, where the values for X max and Y max have to be 
chosen properly. 

The correct configuration space asymptotic boundary conditions for the 3N scattering problem were 
discussed by Merkuriev, Gignoux, and Laverne [24]. It should be obvious that the solution of the 
coordinate space Faddeev equation poses a greater challenge at energies above rather than below the 
three-body break-up threshold. In fact below the threshold one has to consider only the (2+1) channel. 
Then, for the asymptotic solution one may divide the sets {A} in those which have the proper quantum 
numbers to allow for a N-deuteron state and those where this is not possible. The 0l A l of the latter 
have to vanish asymptotically, whereas the 0l A l of the former can be factorized asymptotically into a 
radial deuteron wave function in variable x times a free solution in variable y. Normally the set {A} of 
quantum numbers includes L, S, J, where the channel spin S, with possible values 1/2 and 3/2, results 
from the coupling of deuteron spin jd = 1 and single nucleon spin 1/2. Then S is coupled to L, the 
quantum number for the orbital angular momentum associated with variable y, to total spin J. For 



any state \ J n ) one has either two (J = 1/2) or three (J > 1/2) possible combinations of S and L. 
The asymptotic free solution in y consists of regular and irregular parts, which are described by Bessel 
and Neumann functions, respectively. Their relative admixture is governed by the usual scattering 
parameters, namely by phase shifts and mixing parameters. For J = 1/2 one has two phase shifts and 
one mixing parameter, while for J > 1/2 there are three phase shifts and three mixing parameters. 

The inclusion of a 3NF is rather unproblematic in the Faddeev formalism, since it is a short-range 
force (for details see [22]). More complicated however is the problem of taking into account the long- 
range Coulomb force although how this is accomplished was shown by Merkuriev [25] . 

Distinct from the NN problem, where not more than two partial waves are coupled, in the case of 
the 3N problem one has in principle an infinite number of coupled channels for any \J n ) state. However, 
by restricting the interaction to specific NN partial waves one obtains a finite number of channels. For 
example, if one restricts the interaction only to 1 s and 3 si- 3 di NN partial waves, one obtains five 
coupled channels in a 3N ground-state calculation. The number of channels is increased to 34 if one 
includes the interaction of all NN partial waves with total angular momentum j < 4. 

The generalization of the Faddeev formalism to more than three particles was carried out by 
Yakubovsky [26]. Here we turn directly to the 4N problem in configuration space. In case of four 
nucleons there are two different classes of sets of Jacobi coordinates, namely a K and a H class with 

K K K r i ~r r j ~r v k / 0/1 \ 

*ij = y«,fc = y«,fc , %M = r * 3 ' ( 24 ) 

= x,,, y^ = r I -r fc , ^ = T -h±lL - T -i^Ei . (25) 

Since there are 24 permutations one has 48 different sets. Again the cases with Xy and x^j describe the 
same physical situation, and in addition one has to consider that for the H sets yfc; and yik are also 
equivalent. Thus one remains with 18 (12K+6H) physically non-equivalent Jacobi coordinate sets. 

Similar to the three-body case one splits the wave function \I/ into the FY components ipij. By 
means of the FY components one obtains a set of coupled equations, here 18 equations, which added 
together lead to the four-body Schrodinger equation as detailed below. 

The FY components are defined by 

(xyz|tf) = Ylil>ij(xij,y,z) , (26) 

i<j 

where 

^■(Xtf, y, z ) = V£ (x^ y* fc , z« kl ) + ^ (x^-, y* , z* >fc ) + ^ (x 4J , y=, zg >w ) , (27) 
with k < I. Again, if one requires the six coupled equations (one for each ij pair), 

(T - E)ipij (xjj ,y,z) = -^-J] ipH (Xij , y , z ) , (28) 

k<l 

to be satisfied, one easily verifies that the resulting ^ of Eq. (26) fulfills the 4N Schrodinger equation, 

(T-£ + £^)(xyz|*>=0. (29) 

i<j 

Each of the six equations is further split into three sets of coupled equations. The first set is given by 



(T-E + Vy) V* (xy , yg >Jfc , z« k j) = -Vy ^* (x ife , yf fcj , zf fcjV ) + ^ (xtt, yf fe)i , j) 

+ i'fki^-jk-, yfk,v z fki,i) + i^jki^jk: yfk.ii z fki,i) + ^( x ifc, y^, z^^j) + ^f k {yijk, yfi, zfk,il) 



(30) 
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For the second set one has the same six relations (30) but with k ■<-)■ I, while the third set reads 



(T-E + Vytyy (xy , y* z* H ) = -V l3 (x M , y^, z^) + jfc (x«, y^, z*. i) + yg, 

(31) 

If one of the components ^(xa^J^jZ^J should appear in one of the right hand sides of the 
three sets of equations above with A > (/,, it has to be understood as being identical to that obtained 
by exchanging A and ji. 

Similar to the Faddeev case with three nucleons it is not necessary to solve the full set of coupled 
equations for the case of identical nucleons. In fact using permutations one obtains from one ip K and 
one ip H component the other eleven ip K and five ip H components, respectively. Thus one remains with 
a system of two coupled equations (for further details see [27, 28, 29]). 

In order to solve the remaining two coupled equations, and in analogy to Eq. (23), one makes the 
following expansion of the FY components in channels {A}: 

^/*(x, y, z ) = £ **< KiX ' V ' Z) Fli%(x, y, z) , (32) 

{a} xyz 

where, again, for simplicity the subscripts of variables x, y and z and the dependence of the functions 
%p H / K and Fjff K on spin and isospin variables are dropped. As for the 3N case one has the regularity 
condition for bound and scattering states 

<f>H/K( x = °> f > z ) = <I>h/k( x > V = 0,z) = 4>h} k (x, y, z = 0) = . (33) 
For bound states the asymptotic boundary condition reads 

4%} K {x = X max , y, z) = $} K (x, y = , z) = $} K (x, y, z = z£?) = (34) 

with proper values for X max , Y^P^ , and Z^. For the asymptotic boundary condition the situation 
is more complicated. In case of (3+1) scattering below the first inelastic threshold one has for </>j^ the 
same boundary condition as for the bound-state case, and for one has to divide the sets {A} in 
those where a factorization into a 3N bound-state Faddeev component with total angular momentum 
j xy = 1/2 and a free solution in variable z is possible, while the other sets have to vanish asymptotically. 
As for the 3N case the free solution consists of regular and irregular parts combined with proper factors 
containing scattering parameters. The channel spin S K , which results from the coupling of j xy and 
single nucleon spin 1/2, can be equal to or 1 and couples with the orbital angular momentum quantum 
number L associated with variable z to total spin J. For a given state \J W ) one has either one phase 
shift (J = 0) or two phase shifts and one mixing parameter (J > 0). For the boundary conditions in 
the presence of an open (2+2) channel or further break-up channels we refer to Ref. [28]. 

As was the case for the 3N problem a 3NF and in particular the Coulomb force can also be built 
into the FY r-space formalism for the 4N problem (see [30, 31, 32]). However, unlike the 3N case, a 
restriction of the interaction to specific NN partial waves does not lead to a reduction of the number of 
coupled channels. In an actual calculation it is of course impossible not to make a truncation, which 
then needs to be carefully checked for convergence. In order to better examine the circumstances we 
need to introduce a coupling scheme for the angular momenta of K and H class coordinates. They are 
given by, {([(l«s«)j«l/2]S K l«)j^(l? 1/2)^} J and {[(W)j^f }S H (l»s»)j»}J, where l$ K (u = x, 
y, z) is the angular momentum quantum number associated with variable u, while s^ K denotes the 
total spin of the particle pair with relative coordinate u. 

As has been observed in Ref. [33] one can limit the number of states in a useful manner. For fixed 
1% and fixed J and total isospin T the number of K states is strictly finite once the NN interaction 



o 



1: Number of channels for given values of maximal inter-cluster orbital angular momentum 
and maximal two-body angular momentum = j m ax (from Ref. [33]). 
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Table 2: 4 He binding energy in MeV for various NN potentials (no Coulomb force) from FY calcula- 
tions [1, 34]. Results of columns 2-4 include interaction only in NN channels 3 Si- 3 rf 1 and 1 So, whereas 
for those of column 5 the full NN potential is taken setting l^lUx — 1 and j^b^x = 3 (see also Table 1). 



NN Potential 


only 3 Si- 3 di and 1 sq interaction 


full potential 


Number of Channels 


5K+5H 


15K+9H 


27K+9H 


98K+90H 


Nijmegen [36] 


24.16 


24.55 


24.53 


25.03 


Paris [37] 


23.20 


23.63 


23.60 


24.26 


AV14 [38] 


23.36 


23.77 


23.75 


24.62 



is assumed to act only up to j£ = j max - In the same way the number of H states is strictly finite for 
fixed and j£ < j max . The number of channels for various choices of j max and maximal values l^J^ 
is illustrated in Table 1 for J = and T = 0. It is evident that these numbers grow quite quickly with 
increasing j max and Z^ x . 

In momentum space the calculation proceeds along similar lines as in coordinate space, namely by the 
determination of FY components. The principal difference lies in the fact that the integral equations 
automatically take care of the asymptotic boundary conditions. On the other hand, in momentum 
space calculations, one has to work with t-matrices which contain singularities and which have very 
complex structures (see Ref. [22]). Results for the 4 He binding energy via rigorous FY momentum 
space calculations have first been obtained by Kamada and Glockle [1, 33, 34, 35] in the early 1990s. 
In Table 2 we display a few of these results for what were realistic NN potential models at that time. 
Though the number of channels is infinite even if one restricts the interaction to the NN partial waves 
3 Si- 3 g?i and 1 So ) one sees that for this case 15 K and 9H channels are sufficient to obtain a convergent 
result. Considering the full interaction, however, many more channels are needed as can be implicitly 
deduced from the last column in Table 2, 

In more recent FY bound-state calculations the precision has been further increased and the trun- 
cation scheme changed somewhat. Ref. [39] employed modern realistic charge dependent NN potential 
models (plus additional 3NFs) and the following channel truncation scheme: < 6 and < 8, 

as well as J2u — 14 an d J2 U — 14- The most sophisticated calculation included couplings to T = 1 
and T = 2 isospin channels and required a total of 4200/T and 2000iJ channels. 

For the scattering problem truncation schemes different from that of Table 1 have been used. In 
order to describe n- 3 H and p- 3 Re scattering below the three-body break-up threshold with a precision 
of about 1% it is necessary to include all channels with j x < 4, j y < 4, and j z < 3 [40] in a j-j coupling 
scheme (e.g. j-j coupling for K sets: {[(l^s^(l^l/2)j^]j^ l/2)jf } J). 

As pointed out above the Yakubovsky scheme has been formulated for any number of particles, but 
this scheme has not yet been used in calculations for A > 4 systems. A first more formal application 
for A = 6 has been given in [41], where also some hints for a numerical implementation are provided. 

2.2 The Alt-Grassberger-Sandhas (AGS) Formulation 

Starting from Faddeev and FY theory for treating the many-body scattering problem Alt, Grassberger 
and Sandhas introduced a different set of equations, which today are known as the AGS equations [42]. 
However, apart from just working with a different set of equations, the main motivation of practitioners 
of the AGS formalism was to seek possible further reductions of the problem. In fact the assumption of 
separable forms for the NN t-matrix in the resulting AGS subsystem amplitudes leads to one-variable 
integral equations which are much simpler to solve than the Faddeev or FY equations. 

We describe the AGS formalism by considering the 3N scattering problem for the (2+1) partition 
as introduced in the first part of Section 2.1. The central point of the AGS formalism is the calculation 
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of a transition amplitude defined in analogy to the two-body case 



(*p\Ufr\* a ) = (Q P \VP\9P) , (35) 

where <3> a , V 13 , and are defined as in Section 2.1. Considering first the case U aa one has 

U aa \$ a ) = V a \¥ a ^) = (yp + VJlVP). (36) 

Using Eqs. (6) and (7) for one finds 

C/ Q a|$a> = V£G^|tt«) + V^V^)) = VpG pU p a \<S> a ) + V^U^ a ) . (37) 

In a similar way one obtains for Up a with (3 ^ a: 

Up a \$ a ) = V a \<$> a ) + V a G a U aa \$ a ) + ^ 7 G 7 C/ 7Q |$ Q ) . (38) 

In Eqs. (36-38) a, (3, and 7 are all different from each other and can take the values 1, 2, and 3. One 
should note that for a given entrance channel a is fixed resulting in one equation of type (37) and two 
equations of type (38). Thus, as in the Faddeev case, one has a set of three coupled equations for 
U aa , Up a , and U ia . As for the Faddeev components, one may benefit from the fact that the nucleons 
(neglecting the electromagnetic interaction) can be considered identical particles, and reduce the AGS 
equation set to a single equation. Note that the operators Up a act in the full Hilbert space and carry 
information of the three-body break-up channel as well. 

In order to bring the AGS equation set to its standard form one uses 

Va\<$>a) = G 1 \§a), V a G a = t a G , (39) 

where t a is the NN t-matrix, however, embedded in the 3N-space. Combining Eqs. (37) and (38) into 
a single equation one finds 

U{} a \& a ) = G^Sapl^a) +^2S^t y G U 7a \^ a ) , (40) 

7 

where 5 a p — 1 — S a p. Therefore the standard AGS equation for Up a reads: 

Up a = G Q ^ap + y^Sfryt^GoUyg . (41) 

7 

As in the case of the Faddeev equations an incorporation of a short-range 3NF into the AGS equations 
does not lead to serious complications (see e.g. Ref. [22]). More difficult, however, is the incorporation 
of the long-range Coulomb force in momentum space, although practical high-precision solutions turn 
out to be surprisingly simple [3] (instead of a Yukawa cut-off a smoothed step-function cut-off is taken). 

In order to discuss the AGS equations for the 4N system we first introduce some useful notation. In 
analogy to the 3N case we encode the six possible pairs [(1,2), (1,3), (1,4), (2,3), (2,4), (3,4)] by Greek 
letters. In addition we encode the sub-cluster partitions, which are four of the (3+1) type [(ijk)l, (lij)k, 
(kli)j, (jkl)i] and three of the (2+2) type [(ij)(kl), (ik)(jl), (il)(jk)}, by capital Latin letters. 

The AGS formalism for the 4N system is reviewed in great detail in Ref. [43]. There it is shown 
that the two-cluster sub-amplitudes Up a satisfy the equation 

U L p a = G^ap + E Ifht^GoU^ , (42) 

7 

where pairs a, (3 and 7 must belong to one of the two sub-clusters defined by partition L. Thus, in case 
of a (3+1) partition they must all belong to the three-body sub-cluster. A comparison with Eq. (41) 

1 1 



for the 3N system, shows that Up a has to be understood as the AGS 3N t-matrix operator embedded 
in the 4N space. In the case where L encodes a (2+2) partition Up a is the AGS t-matrix operator for 
the scattering of two noninteracting pairs. 

As in the 3N case, in the 4N case an analogy with two-body scattering theory can be made by 
defining a transition amplitude U^ L such that the matrix elements of U^ L give the appropriate two- 
cluster to two-cluster transition amplitudes for L — y M scattering. As shown by Alt, Grassberger, and 
Sandhas [42] one obtains the following set of equations 

= {GotpGo) 1 5lM^aP + $ MNUfaGpt^GoU^ , (43) 

V7 

where pair a ((3) must belong to one of the two sub-clusters defined by partition L (M). Note that, 
due to the definition of U^, pair (5 must also belong to one of the two sub-clusters defined by partition 
N, which, due to 5mn, is different from partition M. Thus M and N cannot be both (2+2) partitions. 
Therefore one may split Eq. (43) into two equations 

(GotpGo) 1 5KK'8 a (3 + X! ~&K'K"Up K 1 Gq^GqU^ k + ^ U^GqI^GqU^ , (44) 

K"-y H~i 

(G tpG y l 5 a p + ]T U^'G t 1 G U^' K , (45) 

where H represents (2+2) partitions, while K, K' , and K" denote (3+1) partitions. Obviously we 
have set L = K (initial channel represents a (3+1) partition). Of course the setting L = H is also 
possible and leads for an initial (2+2) channel to a similar set of equations. Selecting a specific Ug a K 
one obtains via Eq. (44) couplings to three other transition amplitudes U^' K and to one transition 
amplitudes , which then lead to new couplings via Eqs. (44) and (45). Iterating this process one 
finally arrives at 12 different equations of type (44) and 6 different equations of type (45) which together 
form a set of 18 coupled equations. In analogy to what happens for the FY components described in 
Section 2.1, one has to solve 18 coupled equations for the two-cluster transition operator. 

As for the FY case the 18 coupled equations can be reduced further for identical particles. In 
fact, using permutation operators the equation set is brought down to two coupled equations with two 
transition amplitudes to be determined, one of the U K ' K and one of the U HK type. In case of an initial 
(2+2) partition state one ends up with two coupled equations of types U KH and U H ' H . 

In order to evaluate scattering observables one has to sandwich the transition operators of Eqs. (44) 
and (45) between the respective components of the initial and final state two-cluster wave functions. 
For the (3+1) partition the FY components are constructed from standard Faddeev components of the 
bound cluster, spin and isospin of the single nucleon, and a plane wave for the relative motion. For the 
(2+2) partition the FY components are constructed from the deuteron wave functions. 

Even though the AGS technique aims at a direct determination of the scattering amplitudes the 
formalism may also be used to calculate bound-state wave functions. In this case one has to solve a 
system of 18 homogeneous equations for FY amplitudes, 

W = T, G otaG UipS LM \il!f), (46) 

MP 

which again can be reduced to a system of two equations by considering identical nucleons. 

The AGS equations do not depend on a specific representation, but they have been worked out for 
momentum space calculations only. In this case the relevant equations are projected onto a complete 
set of momentum states, which are expressed using Jacobi coordinates. In the 4N case this leads to 
three momentum variables plus additional spin and isospin variables. Finally, projecting onto a specific 
state \J n ), one obtains an equation system with an infinite number of channels in three continuous 
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Table 3: 4 He binding energy in MeV for NN potential model Y4 [44] obtained in the AGS calculations of 
Ref. [44] for various ranks N in the separable ansatz of the (3+1) subsystem amplitudes (only j w = 1/2+ 
for the 3N subsystem is taken into account). 



rank N 


N=l 


N=2 


N=3 


N=4 




32.370 


32.368 


32.675 


32.666 



variables. As already discussed for the FY case of Section 2.1 one then has to truncate the set of 
equations appropriately while checking the convergence behavior of the desired observable. 

The reason for the preference of using momentum space for AGS calculations is due to the fact that 
powerful approximations can be be incorporated quite easily. As already pointed out in the beginning 
of this section one can employ separable non-local NN potentials. Their use in Eq. (43) leads to new 
equations for so-called subsystem amplitudes, which now depend on two rather than three continuous 
variables. Such an approach is usually labeled as [2V] in AGS calculations. In a [2V] calculation one 
remains with a problem similar to the Faddeev case for three particles interacting via local potential 
models. However, now the number of coupled equations equals eighteen rather than three. One can 
further reduce the problem by making an additional separable ansatz, namely for the subamplitudes, 
which then leads to one continuous variable only (usually labeled as [IV]). The separable ansatz for 
the subsystem amplitudes should in principle be brought to convergence by increasing the number of 
terms, the so-called rank, of the expansion. This however leads to an increase in the number of coupled 
channels to be considered. 

In principle it is not necessary to use separable expansions for the (2+2) subsystems. The set of 
equations (44) and (45) can be reformulated such that the contribution of the (2+2) subsystems is 
calculated exactly. This then leads to the so-called convolution method, which, as stated in Ref. [43] , is 
particularly advantageous when used simultaneously with the [IV] formulation for the (3+1) partitions, 
since in this case one has only to calculate an extra integral over the singularities of a resolvent. Such 
a procedure is usually labeled as a [1V+C] calculation. 

The AGS [1V+C] approach had been considered to be quite promising. In fact the first 4N calcula- 
tion with a full consideration of the tensor force in the 3 s\- 3 di NN partial wave has been carried out in 
this way [44]. Various versions of Yamaguchi-like potentials "Yn" for 1 s and 3 s 1 - 3 di NN partial waves 
have been used, while the 3N subsystem has been restricted to the states j w = 1/2+, 3/2 + . Various 
points have been investigated, among them the convergence concerning the separable ansatz of the 
(3+1) subsystem amplitudes, which we report in Table 3. One sees that a rather good convergence is 
obtained with rank N = 4. A comparison to results of [2V] calculations [45, 46], where the tensor force 
had been taken into account only perturbatively via the so-called too approximation, shows a tensor 
force effect of up to about 0.5 MeV. Later when the same potential models were considered in the above 
mentioned FY calculation [33], only small differences of the order of about 0.2 MeV were observed. 
However these calculations did not employ exactly the same channels. 

Although modern realistic NN potential models are not of separable form one may use a separable 
representation of the NN t-matrix. Usually rank 1 representations are preferred in order to prevent an 
increase of the number of channels. As is discussed in Section 3.2.3 it has been found in the last decade 
that such rank 1 representations are not sufficiently precise for all partial waves of the NN t-matrix. 

2.3 The Variational Formulation 

A large variety of approaches to the many-body quantum mechanical problem make use of the Rayleigh- 
Ritz variational theorem [47, 48]. This theorem is applicable if the solution of a useful equation renders 
stationary some proper functional. In particular, one can show that the solution of the Schrodinger 
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equation (for a normalizable state) renders stationary the Rayleigh quotient for H, i.e. the energy 
functional 

m = (W) (47) 

A lemma that complements the fundamental variational theorem states that the value of the energy 
functional with any trial function is always greater than the ground-state energy and equal to that, 
when the trial function coincides with the exact ground-state wave function. Therefore if one wants to 
calculate the energy of a system the equation to solve is 

5E[V]=0, (48) 

that for the ground-state energy corresponds to a minimization problem. 

Numerous approaches use this variational principle to find the bound-state energy of a many-body 
system. In order to maximize the efficiency of the approach the trial function must have a parametrized 
functional form that is both convenient and suitable to the problem to be solved. With it one calculates 
the energy functional and searches for the parameters that minimize it. 

The different approaches that work with Eq. (48) can be divided in two groups according to the 
criteria followed for the choice of the trial function. To the first group belong those approaches where 
the trial function fy? is written as an expansion on a complete (over-complete) set of square integrable 
functions n that respect the symmetries of the Hamiltonian: 

N 

|*r(AO> = £c n |0 n >. (49) 

n=l 

In this way the minimization procedure corresponds to finding the solution of a (generalized) eigenvalue 
problem 

{H-ES)c = 0, (50) 

where H and S are N x N Hermitian matrices of the Hamiltonian (H nm = (0 n |if|0 m )) and overlap 
integrals of the basis functions (S nm = {</>n\<f>m)) , while c is a ^-component vector of the linear param- 
eters c n . With growing N the size of the Hamiltonian matrix represented on the chosen basis increases 
and the true ground-state energy is approached from above. In principle the true result would be ob- 
tained only for an infinite number of basis functions. The basis can be complete as the hyperspherical 
harmonics (HH) or the harmonic oscillator (HO) basis, or over-complete as e.g. in the method de- 
scribed in Section 2.3.2. In some cases the trial function is constructed on a basis that does not respect 
the symmetry of the Hamiltonian, and the latter is restored after diagonalization. In these cases the 
convergence is considerably slower than in the properly symmetrized cases. In general, however, one 
has to consider the trade-off between the much larger basis and the complications in symmetrizing the 
basis functions (see for example the non-symmetrized HH approach of Section 2.3.1 or the Fermionic 
molecular dynamics technique of Section 2.3.5). 

To the second group belong those approaches where the trial function is chosen according to more 
physical considerations. For example it can be suggested by a cluster picture of the system like in the 
resonating group method (RGM, see Section 2.3.4), or by the form of the potential as in the variational 
Monte Carlo (VMC) technique described in Section 2.5.4. In the stochastic variational method (SVM, 
see Section 2.3.3) the variational procedure does not proceed systematically by the diagonalization of a 
larger and larger Hamiltonian matrix, but in a stochastic way (trial and error). 

It is clear that the approaches belonging to the first group have the advantage of being more 
controlled, in fact the increase of the basis is determined by the growing values of specific quantum 
numbers and in case that convergence for the bound-state energy is reached, one can consider it as a 
precise ab initio result. 



The variational approach also allows one to obtain a wave function that can be used to calculate 
bound-state observables other than the energy. However, one has to remember that the difference 
between the exact value of the energy and that obtained with the trial function which minimizes 
the energy functional is an infinitesimal of higher order than the difference between the true wave 
function and ^t- 

A comment concerning excited states is in order here. While with Faddeev and AGS techniques 
the calculation of an excited bound state follows the same line as that of a ground state, i.e. the 
corresponding equations have to be fulfilled for the energy of the excited state, the situation is different 
for the variational approach. By diagonalizing the Hamiltonian matrix on a given finite basis of square 
integrable functions one gets a spectrum of N eigenstates with eigenvalues E { . If the trial function 
carries the quantum numbers of the ground state of the system, then, the lowest E { is the ground-state 
energy and the ground-state wave function is described by the corresponding eigenstate The other 
solutions correspond to excited bound states of the system that carry the same quantum numbers as the 
ground state. A converged calculation for the energy of such excited states often requires more basis 
states than a ground-state calculation. In most cases, however, excited bound states have quantum 
numbers different from those of the ground state. For a calculation of such states one has to implement 
the proper quantum numbers into the basis set when solving the eigenvalue problem. 

For few-nucleon systems only very few, if any, of the excited states obtained by diagonalization, 
actually correspond to bound states; most of them lie in the continuum (Ei > E t ^). Since the basis set 
$ n consists of square integrable functions the do not have the proper continuum boundary conditions. 
In fact proceeding in such a way one obtains a fake discretization of the continuum. 

A variational principle similar to that for bound states has been developed by Kohn [49] for the 
phase shifts and elements of the scattering matrix in nuclear collisions. Following the original reference 
the principle can easily be illustrated for the one-dimensional problem. Considering the collision in 
s-wave of two particles of mass m interacting by a short-range potential with range Rq, one has to solve 
the following equation 

r d 2 

+ k 2 - mV 12 (r) u(r) = , (51) 



.dr 2 

where k is the wave number of the relative motion and Vu{r) is the interaction potential. The radial 
wave function u(x) fulfills the regularity condition u(0) = 0, while for x ^> Rq one has 

u(r) = A sin(Ax) + B cos(kr) = A[sm(kr) + tan 5 cos(kr)] = (A 2 + B 2 ) 1/2 sin(kr + 5) , (52) 

where A and B are some constants and the phase shift 6 is defined by tan S = B/A. If L = u'(a)/u(a) 
indicates the logarithmic derivative of u(x) at x = a ^> R , then u satisfies the equation 



I 



' du{r, , 2 2 



+ k u (r) — mVi2(x)u (r) 



dr + Lu 2 (a) =0. (53) 



\ dr 
o v 

For the binding energy the basis for the variational principle lies in the following: for a given value 
of L, k 2 , as calculated from Eq. (53), is stationary. The crucial point is that one can equally well 
regard Eq. (53) as providing a stationary expression for L for a given value of k 2 . Consequently, 
this is the variational approach more natural in collision problems, where the energy of the system is 
fixed. Having determined L from its stationarization one obtains the phase shift 8 from the relation 
L = kcot(ka + 5). For the general case of a two-fragment scattering the formulation of the variational 
approach just described for s-wave scattering in one dimension may appear more involved, but the 
essence is the same. (In the literature one also finds different versions of Kohn's original formulation, 
see e.g. Ref. [50].) 

In principle, both bound-state and two-fragment scattering calculations (e.g, elastic scattering of 
n- 3 H or <i- 4 He) could be tackled with any of the various variational approaches discussed in the next 
subsections. However, scattering calculations have been carried out only for some of the approaches. 
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In the following we provide a description of various present-day variational approaches. The discus- 
sion of two of them will be, however, postponed to the section devoted to the Monte Carlo technique 
(see Sections 2.5.3 and 2.5.4). 

2.3.1 The Hyperspherical Harmonics (HH) Method 

The HH method is a variational method where the trial function is written as an expansion on the 
HH basis. The hyper-spherical harmonics are the A-body generalization of the spherical harmonics Y\ m . 
After Borel [51] who first employed hyperspherical coordinates in physical problems, it was Gronwall [52] 
who performed an expansion of the Helium atom wave function in HH (though not referring to them 
as such). The first application in nuclear physics is due to a Ph.D. student of Schwinger, R.E. Clapp 
who calculated in his thesis the 3 H ground state with central and tensor forces [53]. 

The HH basis is intrinsically translational invariant. In fact the functions are expressed in terms of 
the hyperspherical coordinates which are defined by a transformation of the N — (A — 1) Jacobi vectors 
...,rj N ). Unaffected by the transformation are the 2N polar angles 9i and fa of the Jacobi vectors 
(rji = fa)). The remaining N hyperspherical coordinates, the hyper-radius p and (N — 1) hyper-angles 
a n , contain information from at least two Jacobi vectors as can be seen from their definition: 

Pn = E^ 2 ; sina n = 1^1, n = 2,...,N; p = p%. (54) 

i=l Pn 

Expressed in hyperspherical coordinates the A-body kinetic energy operator is a sum of two terms 

T = T p + T K (p), with Tp = -^-A p , T K (p) = ^-^f, (55) 

2m 2m p z 

where T p and Tx(p) are given for a system of identical particles of mass m and using mass weighted Ja- 
cobi coordinates [54] . It is evident that this form, with a p-dependent Laplacian and a /iyper-centrifugal 
barrier, is in perfect analogy to the three-dimensional case. 

The %per-angular momentum operator Kjv, which depends on all the (3N — 1) angles (denoted 
by Cljy), is rather complicated so we refrain from giving an explicit expression here. The role of Kjv 
in an A-particle system is an extension of the role of the relative angular momentum operator li in a 
two-body system. The HH basis set 3 ; [x JV ](^7v) consists of the eigenfunctions of 

K%y [KN] ((l N ) = K N (K N + 3N- 2)y [KN] (n N ) , (56) 

where the eigenvalues are expressed in terms of the quantum number and the subscript [Kn] stands 
for the total of (3iV — 1) quantum numbers; for N = 1 one easily makes the following identifications: 
Kf = If, Ki = Zi, y^ Kl ] = Y hmi . 

Besides the y[K N ](&N) one needs also basis functions for the hyperradial part of the wave function. 
Usually one takes expansions in Laguerre polynomials. 

For a better understanding of the HH and their quantum numbers it is instructive to construct them 
recursively, by starting from the N — 1 case discussed above and adding one particle at a time. Following 
the original paper [55] one realizes that the set of quantum numbers characterizing these functions is 
indeed a set of (3N — 1) quantum numbers, which consists of: (i) the angular momentum quantum 
numbers Ii,12---In relative to the N Jacobi coordinates, (ii) the corresponding third components mj, 
and (iii) the (N ' — 1) quantum numbers 112,113. ..tin associated with the degree of the Jacobi polynomials 
that enter the solutions of Eq. (56). 

A useful method for constructing the HH is the so-called "tree" method proposed by Vilenkin [56]. A 
didactic presentation of this method can be found in [57]. Following the tree method, the set of (3N — I) 
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quantum numbers characterizing these functions consist of: (i) the hyperangular quantum number K N , 

(ii) the total angular momentum quantum number L N and the corresponding third component M N , 

(iii) the angular momentum quantum numbers li, I2, ■■■,In relative to the N Jacobi coordinates, (iv) the 
total angular momentum quantum numbers of the n-body subsystems L^_i, Ljv_ 2 ...-^2 (note that L\ is 
missing because it coincides with li), and (v) the hyperangular quantum numbers K N _i,K N _ 2 ---K2 of 
the n-body subsystems. Note that the following relation holds: Kn = 2un + ifjv-i + ^v- 

Constructing the basis functions in such a recursive way makes the evaluation of any two-body 
operator easy, since only the matrix element of the two-body operator between the last two particles 
needs to be calculated, if the wave functions are properly symmetrized. 

The symmetrization of the HH is one of the two main difficulties of the method. In fact in general 
the HH basis states possess no special symmetry under particle permutation. For the 3N and 4N case 
this is a minor problem since a direct symmetrization of the wave function is still feasible (see e.g. [58]). 
For larger systems, it becomes impractical and there is a need for more sophisticated symmetrization 
methods. One approach to the HH symmetrization problem consists in constructing recursively the HH 
by realizing irreducible representations not only of the orthogonal group 0(3 AT) but also of the group 
O(N) (0(3 A/") D 0(3) Cg) 0(N)). The recursive construction makes it possible to label the HH according 
to the group-subgroup canonical chain 0(N) D 0(N — 1) D ...0(2) [59]. This is expedient, since 0(N) 
is the group of the transformations from one set of Jacobi coordinates to another set, and therefore has 
the A-particle permutation group Sa as subgroup (0(3) <g> 0(N) D 0(3) <S> Sa)- 

The classification of the HH with respect to appropriate group-subgroups allows a model space 
identified by K N to be split into smaller subspaces [60]. An efficient technique has been developed by 
Barnea [57] to solve the reduction problem from 0(N) to Sa- It is the application of this technique 
that has for the first time allowed a HH ab initio calculation for more than four fermions [61]. 

Another interesting non recursive procedure that automatically guarantees the antisymmetrization 
of the many-body wave function has been suggested in [62]. There the HH are expanded on the Slater 
determinant basis of the HO shell model (SM). The expansion coefficients can be found by diagonalizing 
the angular part of the multidimensional Laplacian presented in individual nucleon coordinates [63]. 
The applicability of the proposed method has been in principle demonstrated for the case of the 3 ~ 7 H 
and 4_10 He isotopes [62] and for bosonic systems [63] by using simple central interactions. However 
further improvements are necessary to obtain convergence. 

More recently an approach has been proposed, where the basis is not symmetrized [64] . The physical 
states are identified by studying the symmetries present in the eigenvectors obtained from diagonaliza- 
tion of the Hamiltonian matrix. In principle the drawback of such an approach is the large degeneration 
of the basis. However, this problem is circumvented by expressing the Hamiltonian matrix as a sum 
of products of sparse matrices. This particular representation of the Hamiltonian is well suited for a 
numerical iterative diagonalization, where only the action of the matrix on a vector is needed. Up to 
now this method has been applied to six nucleons with a simple central potential. 

The second difficulty in the application of the HH approach is common to all diagonalization methods 
and concerns convergence. In fact for potentials with strong short-range repulsion the convergence rate 
of the HH basis turns out to be very slow. Two methods have been suggested in order to improve 
the convergence rate. One of them couples the HH expansion with the construction of an effective 
interaction (EIHH) and will be explained in Section 2.4.2. The other one is the so-called correlated 
hyperspherical harmonics (CHH) method first used by Fenin and Efros [65]. Correlation functions /(r^-) 
are introduced that correlate the HH basis according to the characteristics of the potential, 

yiK](Ci N ) n/(^)3W^) , (57) 

where r\j is the relative distance of particles i and j. The function f(r) can be obtained from the 
solution of the two-body Schrodinger equation, since at short distances the role of other particles is 
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Table 4: 4 He binding energy in MeV as function of the maximal value of the grand- angular quantum 
number K from a CHH calculation with the AV14 NN potential (with Coulomb force) [71]. 



K 

-''■max 


2 


4 


6 


8 




20.60 


23.12 


23.71 


23.85 



rather unimportant. Though such correlations lead to a considerable improvement [65] one still needs 
in general a rather large number of HH terms in order to reach convergence. The convergence can be 
improved considerably if one introduces state dependent (see Ref. [66]) or longer range correlations [67]. 
However, the two-body long-range correlations are less under control, because correlations among more 
particles become important. An interesting HH reorganization was made in Ref. [68], where HH have 
been constructed that allow a better treatment of long-range correlations typical for halo nuclei. In 
a toy model of a loosely bound 5 He nucleus the clusterization aspect of the system dynamics could 
be described. Unfortunately this approach has not been applied to real nuclear systems. However, an 
extension of this idea to nuclei with two loosely bound valence nucleons, and in particular to Borromean 
nuclei, would probably make it possible to achieve a proper three-body description of such systems at 
large distances. 

The insertion of correlations complicates the calculation of the Hamiltonian matrix elements. In 
fact the orthogonality of the HH is destroyed and one has to solve a generalized eigenvalue problem 
involving the metric matrix. This implies fully 3A — 3-dimensional integrals, which can be reduced to 
3 A — 6-dimensional integrals using the technique described in [69]. Therefore it can be more convenient 
to use the pair-correlation ansatz with an expansion of the wave function in Faddeev components 

tffai, •••> Vn) = E firMjiV^ •••> Vn) , (58) 

i<j 

(where we have dropped again spin and isospin degrees of freedom). Such a Faddeev approach leads 
to at most 12-dimensional integrals. A calculation along this line has been performed in Ref. [61], 
extending the integro-differential equation approach (IDEA) [70] that limits the expansion to HH of 
lowest order. 

The first HH calculations for the 4 He ground state with realistic potential models and with sufficiently 
convergent results have been carried out by Viviani, Kievsky and Rosati [71]. In Table 4 we show results 
for the binding energy with the AV14 NN potential [38] from their first calculation. Inspecting the table 
one sees that the result is not yet fully converged. In following studies it turned out to be more efficient 
to work with HH instead with CHH expansions. In Ref. [72] an improved result of 24.18 MeV has 
been found in a calculation with about 2500 uncorrelated HH channels. This contrasts the about 100 
channels, which have been employed for the CHH result with K max = 8 of Table 4. In fact for CHH 
calculations much less channels are necessary, while on the other hand one has to cope with the loss 
of the orthogonality of the basis functions. Therefore it cannot be said in general which of the two 
methods is advantageous for a specific problem. We should mention that a fully converged HH result 
of 24.23 MeV for the 4 He binding energy with the AV14 potential is given in Ref. [73]. 

It is an important fact that the HH expansion is not restricted to local potential models, but can 
also be applied to non-local forces [74] . 

2.3.2 The Gaussian Expansion Method (GEM) 

The GEM is a variational approach first introduced to study muonic molecules in 1988 [75]. Shortly 
after, it was applied to calculate the 3N and 4N nuclear binding energies for a simple potential [76]. 
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The trial function is expanded on an over-complete basis. In order to account for the components 
generated by the various possible arrangements of the particles that form the system under investigation, 
one allows all possible choices of different Jacobi coordinates. In this connection one has to distinguish 
between the cases of non-identical and identical particles. In the latter case there is no need to employ 
the Jacobi sets that may be obtained from each other via permutations of particles. Therefore it is 
sufficient to use one single set for A = 3 and for A = 4 two sets, for example those which were described 
in Section 2.1. 

The GEM wave function for a nuclear state | J w ) has the following form [75]: 

\J n }= E 4 XN ^A\^), (59) 

7 ,{A7V} 

where the antisymmetrizer is denoted by A, 7 enumerates the different Jacobi sets, {A} represents the 
angular, spin and isospin quantum numbers of the basis functions |<M AA ^), and the set {N} encodes 
information about the radial basis functions (see below). For a better illustration of the ansatz we 
specialize to the wave function | J 71 ") of the 4N system: 

\r) = £ c r } a + e c { h xn} a , m 

{XN} {XN} 

where K and H stand for the above mentioned Jacobi coordinate sets. The basis functions are described 
for the K and H sets by 

(Xf 2 ,yf 2 ,3,zf 23 ,4l$r } ) = t&Xx&t&K^ (61) 
(xf 2 ,y3Vf 2 ,34l*r } > = «W 2 ) &Vu) C«34) 4^2,^4,^,34), (62) 

where we have suppressed the spin and isospin degrees of freedom in the functions F^J K . The Jacobi 
coordinates used in the two equations above are defined in Eqs. (24) and (25) and the quantum number 
set {A} indicates a specific channel (due to the different Jacobi coordinates the quantum numbers for 
K and H sets are different). The radial functions are given by 

<f,W(r) = r l e^ r/r - ) \ (63) 

where / belongs to the set of quantum numbers {A}. Furthermore we have the following equality 
{iV} = {rii, n 2 , n 3 }. The Gaussian range parameters are chosen to lie in a geometrical progression 
(j'n = Tia n_1 , n = l...n max ). Choosing the range parameters in this way allows one to describe both the 
short-range correlations and the long-range asymptotic behavior precisely. The use of Gaussians allows 
the analytical calculation of the Hamiltonian matrix elements. 

In the GEM approach one proceeds systematically enlarging the basis by increasing the values of 
the different quantum numbers until satisfactory convergence is obtained. In other words the number 
of channels A has to be increased to attain sufficient convergence. Of course also the various rii have to 
be chosen sufficiently large. 

One notes that the GEM calculation exhibits some similarities to the FY bound-state calculation 
described in Section 2.1. In fact, H- and i^-type Jacobi coordinates and the various quantum number 
sets {A} are actually identical, the difference being in the calculation of the radial parts of the wave 
function. In the GEM case the radial parts are obtained from a variational calculation, where the basis 
consists in a product ansatz of radial functions of single Jacobi coordinates, as opposed to the FY 
approach, where the radial parts are determined by a numerical solution of the FY equations. 
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2.3.3 The Stochastic Variational Method (SVM) 



The SVM is an approach that was first introduced in 1977 by Kukulin and Krasnopol'sk to calculate 
the three- and four-body nuclear binding energies for simple potential models [77] . About twenty years 
later the first calculations with realistic potentials were carried out by Varga and Suzuki [78]. 

The strategy of this variational method is somewhat different from the GEM method discussed in 
Section 2.3.2, although the trial function is here also expanded in Gaussian functions. For the SVM 
approach correlated Gaussians are used. That is, the arguments of the Gaussians are not simply in terms 
of single Jacobi vectors, but in addition products of different Jacobi vectors. Particular emphasis is 
given to the correlation that must be flexible enough to describe the full variety of correlations between 
the nucleons, e.g., the short-range correlation due to the strong repulsive force, the clustering typical 
in some light nuclei, or the long-range correlation in light halo nuclei. Correlated Gaussians imply 
in general many non-linear parameters. Then again the variational basis is nonorthogonal and over- 
complete, i.e. none of the components is indispensable, and one can replace a component by a linear 
combination of others. This fact actually gives an excellent opportunity to use a stochastic optimization 
procedure. In fact, for a large Gaussian basis the hypersurface of the energy function is generally too 
complicated to locate an absolute minimum of the energy with respect to the Gaussian parameters. 
However, it is usually enough to find a sufficiently low point in energy. In an actual calculation one has 
to decide whether one wants to spend much computational time in a search for the global minimum of 
a smaller basis set or, as an alternative, to add more basis functions in order to lower the energy this 
way. Finding the right balance between both choices is quite a difficult task. 

In the SVM the wave function of the A-nucleon system \J n ) expanded in correlated Gaussians reads 
as follows 

(^,...,^1^) = £ c^A{e^< X]x Fi X \^ fj^)}, (64) 

{A}fc 

where spin and isospin degrees of freedom are dropped, x stands for the set of (A — 1) Jacobi vectors 
r] i , while {A} encodes the quantum numbers of the various channels. In practice the index n in Eq. (49) 
corresponds here to the set ({A}, k). 

The matrix B k has the dimension (A — 1) x (A — 1) and is a positive definite symmetric matrix of 
non-linear parameters B k = B\y The explicit form of the multiplication of vectors x and matrix B k is 
defined by 

x*B k x = f) B% 77, • Vj • (65) 

ij 

Note that the trial wave function is not written as a separate product of angular and radial parts. 
Even though there is an expansion on angular momenta (contained in the set {A}), the angular part 
of the trial function is not exhausted by this expansion, since there are also angular correlations in the 
Gaussians. This means that the angular basis is not a complete basis, but an over-complete one. 

As already pointed out it is not a simple task to optimize this large number of parameters towards 
a minimum. Fortunately here one works with an over-complete basis and thus one may use a stochastic 
approach. The stochastic procedure of the SVM is organized in three steps: (i) one generates a trial 
function by choosing the parameters B^ and the channels randomly; (ii) one judges its utility by the 
energy gained by including it in the basis, and either keeps or discards it; (iii) one repeats this "trial 
and error" procedure until the basis set leads to convergence. 

For a better search of the minimum after a certain number of new basis functions are accepted, one 
may make a cyclic optimization, where one optimizes the entire basis set by tuning the parameters of 
only one function at a time [79]. To optimize the entire basis set in such a way is numerically much 
cheaper than a simultaneous optimization of all the Gaussian parameters of all basis functions. 

In Fig. 1, a result reported from Ref. [78] illustrates an example of the convergence rate. The 6 Li 
binding energy with the Volkov potential [80] has been calculated several times starting from different 
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Figure 1: Results of a SVM calculation [78] for the convergence pattern of the 6 Li ground-state energy 
with the Volkov potential [80] on two different random paths. 



random points and two of these calculations are shown in the figure. The energies obtained from different 
random paths approach each other after a few initial steps and converge to the final energy. The energy 
difference between two random paths as well as the tangent of the curves gives some information on 
the accuracy of the method with a given size of the basis. Besides the large number of non-linear 
parameters, the treatment of the increasing number of channels in the expansion of the wave function 
also poses a formidable task. Therefore one can use an alternative approach to cope with this problem. 
Instead of using the partial wave expansion, the angular dependence of the wave function is represented 
by a single solid spherical harmonics whose argument contains additional variational parameters 

G LMl (x) = v 2K+l Y LMl (v) (66) 

with v = J^i' 1 UiXi. Only the total orbital angular momentum appears in this expression and the 
parameters u,i (as well as K) may be considered as variational parameters. The factor v 2K+L plays 
an important role in improving the short-range behavior of the wave function. This form makes the 
calculation of the matrix elements for nonzero orbital angular momentum much simpler. 

A substantial improvement in the optimization procedure of the stochastic variational approach has 
been made by the analytic gradient optimization method (AGOM) [81, 82]. In order to speed up the 
search of the minimum, in the AGOM one calculates the energy gradient with respect to the Gaussian 
parameters in an analytical way. By indicating with a t one of the many non-linear parameters on which 
the basis function <p t depends, and by elaborating on the differential of the secular equation (50), in 
Ref. [82] it is shown that the derivative of the total energy with respect a t is given by 
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Therefore, by calculating such a derivative for each a t one can get an analytic expression for the whole 
energy gradient. Further details on the analytic gradient calculation and on its very efficient numerical 
implementation are laid out clearly in Refs. [81, 82]. 

Up to the present the AGOM has only been applied to atomic and molecular few-body systems. To 
illustrate its quality Table 5 shows results from applying AGOM to the positronium molecule (molecule 
of two positronium atoms). Also shown are results from SVM calculations where the above mentioned 
cyclic optimization has been used. In the AGOM calculation the possibility to determine the entire 
gradient analytically allows the optimization of all basis functions simultaneously. Inspecting Table 5 



Table 5: Convergence of the total energy (in a.u.) of the positronium molecule in the state with L = 1 
and negative parity 
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200 
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300 
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400 


-0.334 408 266 3 


-0.334 407 971 


500 


-0.334 408 295 5 




800 




-0.334 408 177 


1200 




-0.334 408 234 


1600 




-0.334 408 265 8 



one finds only a very small difference between AGOM and SVM results for 100 basis states, whereas for 
200 basis states the AGOM energy is considerably lower than the SVM energy. It is also interesting to 
see that the AGOM result with 500 states is somewhat lower than the SVM energy with 1600 terms. 
Our example shows that the AGOM is a very powerful technique. 

As already said, the AGOM has not yet been applied to nuclear few-body systems. While in atomic 
and molecular systems one has to deal with the simple Coulomb force, the complicated structure of the 
nuclear force could make the calculation much more involved. However, the method is quite intriguing 
and deserves being tested in nuclear few-body problems. 



2.3.4 The Resonating Group Method (RGM) 

The resonating group method consists of an expansion of the nuclear wave function in cluster wave 
functions. In principle these clusters are chosen, based on energy considerations, according to all 
possible fragmentation channels of the compound nucleus [84, 85]. For example, 4 He can be fragmented 
in five possible ways: (j9+ 3 H), (n+ 3 He), (d + d), (d + n + p), (p + p + n + n). In addition one could also 
consider other clusters like (d + d), where d represents the virtual quasi-bound state of the s-wave np 
pair with isospin T=l. 

A typical RGM trial wave function can therefore be written as 
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where Xj identifies the fragment and $^(Xj) describes its internal wave function, / is the number of 
clusters for a given fragmentation and i enumerates the various /-fragment sets (in the 4 He example 
above / would assume only the values 2, 3, and 4). The functions Ff\qJ]j) are the wave functions for the 

relative motion of the fragments and the variables 77^] represent the corresponding Jacobi coordinates. 

It is evident that the maximal possible value for I is equal to A. In few-nucleon physics one normally 
is limited to / = 2 thus considering a fragmentation of the system into sets of two clusters Aj and Bj. 
The consideration of / > 2 is strictly necessary only for scattering states with open /-fragment channels. 

The RGM approach is founded on the idea that the Hamiltonian can always be rearranged as a sum 
of well-defined cluster Hamiltonians 

H X = T X + J2 V jk, (70) 
i<fcex 
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plus the potential terms among nucleons belonging to different clusters as well as the relative kinetic 
energies among the clusters. In Eq. (70) T x is the total kinetic energy of the nucleons in cluster X 
minus the kinetic energy of the CM of the cluster, while Ax denotes the number of nucleons in the 
cluster. Besides the NN potential Vjk, a 3NF can also be included in Hx- Apart from the so-called cluster 
model, where one uses effective cluster-cluster interactions and where the clusters are either structureless 
particles or described by simple models, in the RGM approach one retains the microscopic interaction 
between the constituents of different clusters, and the cluster wave functions, <3>x(l,2, Ax) = $(X), 
are obtained from the same microscopic interaction. 

The RGM strategy consists in a calculation of the various cluster wave functions $(X) as a first 
step. In a second step, the relative motion wave functions are determined variationally. Let us consider 
for example, the simple case of a single fragmentation of the nuclear system into clusters A=Ai and 
B=Bi with trial wave function 

— *4.[$(A) <&(B)F(r])] . (71) 
The factorized form of corresponds to the division of the Hamiltonian into 

H = H A + H B + H' (72) 

with the two cluster Hamiltonians Ha and Hb and the Hamiltonian H' for the relative motion of the 
two clusters. The latter has the form 

H> = T rcl + £ £ Vy , (73) 

where T re i is the relative motion kinetic energy. A 3NF can also be included in Hb and H' . 

The solution of the variational problem for the RGM equations can be obtained in matrix form (see 
e.g. [86]) or by using integral kernels (see e.g. [85]). 

The philosophy of RGM calculations has to be understood as follows. A bound-state wave function 
or the interior part of a scattering state is expanded on a over-complete set of basis functions. As 
shown above such an over-complete basis consists of two parts: the cluster and the relative motion wave 
functions. As a starting point one could consider all possible two-cluster configurations. Additional 
basis states could include virtual excitations of the single clusters since these should have the least 
overlap with the original basis. Let us consider, for instance, a two-cluster state with an arbitrary 
cluster A and a cluster B consisting of a pn pair with isospin T=0. The lowest energy state of the 
pn pair is the deuteron bound state. All virtual excitations of the pn cluster lie in the continuum and 
are characterized by relative energies E-q > and quantum numbers for relative angular momentum I, 
spin s, and total angular momentum j. The virtual two-body excitation simulates a free (1+1) cluster, 
where, however, none of the two nucleons can escape the nucleus. For any set of allowed quantum 
numbers /, s, and j and any Eb > one has in principle a well-defined continuum solution. However, 
since here one is concerned with short lifetime virtual excitations one may use an expansion in a set of m 
localized functions. In fact, in RGM calculations not only cluster bound states, but also cluster virtual 
excitations are determined in this way. Such expansions lead to a discretization of the continuum. 
Thus in our example the pn continuum is represented by m discrete states which can be ordered with 
respect to growing E B . It should be noted that an increase of m does not only lead to an increase of the 
number of states in the nuclear spectrum, but in general also to a shift of the n-th state in the spectrum. 
Having fixed the model space of the various clusters one could use as variational parameter the total 
energy E tot and allow all two-body cluster states $(Aj)$(Bj) with E Ai + Eb % < E tot . Then E tot should 
be further increased until convergence is reached for the observable under investigation. However, in 
order to reach convergence, it might be preferable to allow for certain clusters to have higher internal 
excitation energies than others, although this makes a controlled expansion a bit difficult. 



Table 6: 4 He binding energy in MeV for various model spaces from RGM calculation [87] with the 
BonnR potential [89] (with Coulomb force). See text. Also given is the corresponding GFMC result [90]. 
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^From the discussion above it is evident that in most cases a calculation with the simple ansatz (71) 
will not lead to very realistic results. Other fragmentations should be taken into account in addition. 
To better illustrate the RGM procedure we consider an example from the literature. In Ref. [87] the 4 He 
compound system is studied by employing a realistic Vnn (f° r a more recent calculation with additional 
consideration of a 3NF see Ref. [88] by the same authors). The variational calculation is carried out 
using Gaussians as basis functions. In Table 6 we display results for the 4 He binding energy for various 
model spaces. By including only the two-cluster states p- 3 H, n- 3 He, and d-d, one underestimates the 
corresponding GFMC precise ab initio result for the binding energy by 12 MeV. By adding the d-d 
channels 5 do and the d-d channel 1 so- ls o one gains almost 2 MeV. In order to reach a high-quality result 
additional (i) 30 p- 3 H and 30 n- 3 He virtual excitations, (ii) 82 np-np virtual excitations (both pairs 
with T=0), and (iii) 80 (NN)-(NN) virtual excitations (both pairs with T=l) are necessary. 

It is interesting to note that because of the over-complete basis it is not necessary to determine the 
various cluster wave functions with a very high precision, e.g., in Ref. [87] the deuteron wave function 
is described by only three Gaussians for the s-wave and two Gaussians for the rf-wave. This leads to 
an underestimation of the deuteron binding energy of about 0.3 MeV. Probably more accurate cluster 
wave functions would lead to a reduction of the number of virtual excitations, but, of course, one has 
to search for the best compromise, which makes a convergent expansion of a RGM calculation a bit 
tricky. On the contrary, for the case of a scattering calculation with two asymptotic fragments A and 
B it is desirable to calculate the corresponding wave functions of clusters A and B with high precision. 

2.3.5 Fermionic Molecular Dynamics (FMD) 

The FMD approach has been suggested by Feldmeier [91] in order to solve the many-body problem of 
interacting identical fermions with spin 1/2. The aim has been to describe nuclear ground states and 
heavy-ion reactions in the energy regime below particle production threshold (for a review see [92]). 

In FMD the nuclear wave function is a linear combination of Slater determinants of single-particle 
states \qi) given by linear combinations of Gaussians with complex parameters b k and a k , 

(x|ft> = £c fc exp j- (^r) } lXk) l&) ' (74) 

and additional spin and isospin wave functions \xk) and |^), respectively. For \xk) all orientations in 
spin space are possible (\xk) — xi\ T) + Xk\ D)- While in most of the FMD applications represents 
either a proton or a neutron state, in more recent studies isospin mixing has been allowed [93] i.e. 

i~ I 

16c) — Ck\ t) D- Such a superposition introduces a charge mixing in single-particle and many-body 
spaces. Therefore, in order to conserve charge one has to project onto a state with a fixed value of the 
third component of the isospin for the nuclear system under consideration. The long-range correlations 
in nuclear wave functions originate from a pion exchange between two nucleons. Besides charge transfer, 
due to its pseudoscalar nature the pion induces a parity change. The FMD states can be chosen not to 
be good parity states. Parity is then restored by a proper parity projection operator. 

The energy minimization is made with respect to all the various single-particle parameters. In 
general, the resulting correlated many-body state breaks translational and rotational invariance. Hence, 
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after minimization, the wave function is projected onto zero total momentum and total spin J [94]. The 
FMD formalism in its present form relies on an operator representation of Vnn in coordinate space, 
which limits the number of potential models that can be used. A possible choice would be the AV18 
potential [96]. It has however a rather strong short-range repulsion, which makes it very difficult to 
reach convergence in a FMD expansion. Thus softer versions of AV18 are often used by employing 
unitary transformations with the UCOM method [12]. (Such a potential, however, does not seem to be 
ideal for ab initio calculations as discussed in Ref. [95]). Despite all this, the FMD technique might be 
quite promising and is in principle systematically improvable. Up to the present, however, it is difficult 
to judge the precision of the FMD, since converged benchmark tests are missing in A > 4 systems. 
Thus, here we can only report about a case, where a 4 He FMD calculation has been made with different 
UCOM versions of AV18 taking just a single Slater determinant [93]. In comparison to precise NCSM 
and EIHH ab initio results an underbinding between 12.8 and 3.6 MeV is observed [93]. 

2.4 The Similarity Transformation Formulation 

Another way to reformulate the quantum mechanical many-body problem is achieved with the help of 
similarity transformations [97, 98, 99, 100]. To this end we consider the following mean value 

/ = <*|F|tt), (75) 

where is an eigenstate of the Hamiltonian H. The mean value / is invariant under similarity 
transformations e s , i.e. 

f = (^\ e - s e s Fe~ s e s \^} = (^1^) (76) 

with 

|$) = e 5 |^), |$) = e - 5t |^), J= = e s Fe~ s . (77) 
One may consider a subspace P of the Hilbert space with eigenprojector P defined by 

N 

P=£|0„>(0n|, (78) 
n=l 

where the \<p n ) are the eigenfunctions of some convenient Hamiltonian H . The residual space Q has a 
corresponding eigenprojector Q — I — P. Therefore one can write Eq. (76) as 

/ = (<5|(P + Q)P(P + Q)|$) = ($|PPP + PTQ + QPP + QPQ|$> . (79) 

Given a |<3>) entirely contained in the P-space one has 

/ = ($|PPP|$> (80) 

if the following decoupling condition is satisfied 

QJP = Q e s F e~ s P = . (81) 

This means that by solving the decoupling equation (81) it is possible, in principle, to determine S and 
therefore calculate / as the mean value of the effective operator T on the P-space. 

Since the similarity transformation e s represents an infinite sum of operators it can be problematic 
to deal with it in actual calculations. On the other hand there exists a relation that is very useful in 
applications (see e.g. Section 2.4.3), namely the Baker-Campbell-Hausdorff (BCH) expansion 

T = F + [P, S] + i[[P, S], S] + i[[[P, S],S], S] + i[[[[P, S], S], S], S] + ■ ■ ■ . (82) 



Similarity transformations have been used extensively to find the ground and low-lying state energies 
of a system. In this case one has to solve Eq. (81) replacing F by the Hamiltonian H. Unfortunately 
the full solution of Eq. (81) is a very difficult problem. Various methods have been devised to approach 
the solution by limiting the operator S, which is in general a sum of n-body operators S^ n \ i.e. 

S = S [1] + S [2] + S [3] + ■■■ + S [A] , (83) 

to the first two or three terms (see Sections 2.4.1 and 2.4.3). 

One strategy to satisfy the decoupling condition (81) was given in the work of Okubo [97] and later 
revived by Suzuki and Lee [100]. It consists in requiring the operator S to have the property 

S = QSP , (84) 

which implies PSP = QSQ = PSQ = 0. If S has these properties it allows one to write 

e s = l + S. (85) 

and the decoupling condition (81) is also satisfied. Finding S via the decoupling condition (84) would 
permit one to get ground-state and some low-lying state energies as the eigenvalues of the effective 
Hamiltonian H = e s H e~ s in a convenient P-space. 

In order to work with a Hermitian effective Hamiltonian, the similarity transformation e s must turn 
into a unitary transformation U = e G with G = —G^ = arctanh(S — S^) [97, 101, 102], where again 
S has to satisfy the decoupling condition (84). The authors of this formulation called this approach 
unitary model operator approach (UMOA), which today, however, has a more specific meaning (see 
below). Today the approach is more often referred to as Suzuki-Lee (SL) approach. 

The operator S in Eq. (84) can in principle be determined as follows: (i) let the P space of dimension 
dp be spanned by states \oi p ), (ii) take a set of the lowest d p eigenstates \n) of the Hamiltonian H, such 
that P\n) are linearly independent, (iii) take a set of d q states \a q ) belonging to the space Q, then the 
following equality holds 

dp 

(a q \n) = J2(a q \S\a p ) (a p \n) . (86) 
P =i 

By defining matrices Aq and Ap, with elements (a q \n) and (a p \n) respectively, one finds the matrix 
elements (a q \S\a p ) of S via 

S = A Q A P 1 . (87) 

Since a solution of Eq. (86) at the A-body level is not feasible for larger A, in practical applications it has 
first been solved at the two/three-body level, which then leads to a two/ three-body effective interaction, 
which is used to accelerate the convergence of the calculation. As will be discussed in the following 
this has been done in the NCSM as well as in what is today known as the UMOA approach. In both 
approaches the unitary transformation has been applied to the two/three-body HO Hamiltonian with 
a bare interaction. The difference between the two approaches lies in their procedures for obtaining the 
two-body (three-body) effective interaction, e.g. in the UMOA case the two-body effective interaction 
depends on the HO quantum numbers of the center of mass of the pair. This state dependence makes 
the calculation more computationally intense. Moreover the approach is supplemented by a second 
step in which the effective interaction undergoes another unitary transformation, such that the matrix 
elements of the new effective interaction between two-particle-two-hole excitations reduce to zero. In 
this way the diagonalization of the A-body Hamiltonian with the new effective interaction is made only 
between one-particle-one-hole states. Calculations of heavier systems, like 40 Ca and 56 Ni, have become 
possible in this way [108], even though at the price of a huge numerical effort. 

In the EIHH approach the unitary transformation has been applied at the two- /three-body level. 
Because of the HH basis the effective interaction turns out to have a considerable dependence on the 



A-body system. The complications of the antisymmetrization problem for HH functions, however, has 
limited its applications to A < 8. 

Here we would like to mention another approach that makes use of similarity transformations namely, 
the similarity renormalization group approach. The latter draws inspiration from the ideas of renor- 
malization group theory, which can be paraphrased as follows: when a problem is too hard to be solved 
in one step one can still arrive at the solution by breaking it into small steps. So, instead of solving the 
decoupling condition of Eq. (81) directly in order to obtain an effective Hamiltonian H, one can try to 
build H step by step. 

In the approaches described above % is obtained from a unitary transformation U = e G with an 
antihermitian G. The crucial point of the SRG approach is the introduction of a continuous flow 
parameter a such that 

U(a) = U\a)HU{a), (88) 

however, only when a flows towards a specific value a does one have the proper U which leads to the 
diagonal form of H. The idea is to reach that stage by letting a flow. In fact, U(a) can in principle be 
obtained through solving the differential equation, 

^^=[ V (a),n(a)], (89) 
which itself is obtained by differentiating Eq. (88) with 77(a) being an antihermitian operator: 

v ( a ) = ^Epitf(a) = -77 f (a) . (90) 

There are an infinite number of generators rji, since one can construct an infinite number of different 
Ui(a) that tend to the proper U for some value of a^. However, an interesting observation is that the 
relation 

77(a) = [H d (a),H{a)\, (91) 

where T-Ld( a ) is the diagonal part of the Hamiltonian in a given basis, ensures that the off-diagonal 
elements of H(a) decay exponentially with a [103]. This means that one can solve Eq. (89) for H with 
77 given by Eq. (91). Then, by increasing a step by step in the numerical solution of Eq. (89) one can 
make %{a) more and more (block)-diagonal. 

An important remark is in order here. The generator 77(a) is a many-body operator, therefore it 
will generate an A-body effective Hamiltonian and an A-body effective potential. In most cases it is 
inconvenient or impossible to implement such an operator in a many-body code. What is then often 
done (see Ref. [11, 12] and references therein) is to express Eq. (91) in the two- or three-body space, 

7/[ 2 ]([ 3 ])(a) = [^ 2l([3l) (a),^[ 2 ]([ 3 ])(a)], (92) 

and let H^ 3 ^(a) flow versus a more diagonal form. The resulting two- or three-body operator can 
then be implemented in the many-body code. 

In applications to nuclear forces [11, 12] one normally uses instead of the running T-Ld( a ) a different 
diagonal Hamiltonian matrix, like e.g. the kinetic energy which is diagonal in momentum space. 

Contrary to NCSM and EIHH approaches in Sections 2.4.1 and 2.4.2, where one in principle recovers 
the bare interaction from the effective interaction in the limit of an infinite P^-space, in the SRG 
approach the potential is obtained by truncating the flow at a two- or three-particle level thus violating 
the unitarity of operator U even for an infinite Pa space. Therefore the violation cannot be compensated 
in the subsequent many-body calculation by increasing the R4 space. 

We will end this chapter by mentioning UCOM, a method that allows one to obtain a two-body 
effective potential, phase equivalent to a realistic Vnn, but much softer at short range, and therefore 
suitable to use in A-body calculations. In a way this method is similar in spirit to the SRG method 
when Eq. (89) is solved with 7/ 21 of Eq. (92). 
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2.4.1 The No Core Shell Model (NCSM) 

Two methods have facilitated increasing the domain of the few-body field to nuclear systems with 
A > 4, while working with realistic potentials. These are the GFMC and the NCSM approaches. 

The name NCSM is due to the fact that all the nucleons of the nucleus are active and explicitly 
taken into account as degrees of freedom and thus it is not assumed that there is an inert core. The 
NCSM couples the traditional SM advantages of working in a HO basis with the rigor of an ab initio 
approach. For not so soft potentials convergence can be achieved by using an effective interaction. In 
the 1990s NCSM calculations were made with a G-matrix approach (see e.g. Ref. [104]), whereas in 
Ref. [105] a different and more successful direction was taken, which will be described in the following 
and which are based on the SL transformation formulation explained in Section 2.4. 

There are two versions of the NCSM, which differ in their treatment of the CM degrees of freedom. 
In the first version, the Hamiltonian of Eq. (1) is modified [106] by adding a harmonic oscillator CM 
Hamiltonian Hqm to the intrinsic Hamiltonian H 
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where Rcm is the CM coordinate, while Vj k is a renormalized potential depending both on the HO 
frequency Q and A: 

mVt 2 . , 9 

y 



2A 



iY 



(95) 



Naturally, the added CM HO term has no influence on the internal motion. This means that once 



the ground-state energy E$ of H^ 1 is found, then the ground-state energy E of H is obtained by 
subtraction of the CM ground-state energy 3M1/2. For excited states one has to avoid CM excitations. 
This can be achieved by the replacement Hcm — ^ ^^cm w ^h sufficiently large A (Lawson term). A 
good check on the convergence of E is through the independence of the result on the HO frequency f2. 
On the other hand, one may search for the frequency Q which exhibits the best convergence pattern. 

The calculation of Eq^ is performed in a finite model space P^. This is the space spanned by all the 
A-body HO Slater determinants formed by filling the single-particle HO eigenstates with N < iV max , 
where N is the total number of single-particle HO quanta. The residual space Q^, together with P^, 
exhausts the full Hilbert space. 

The approach that is chosen to construct the effective Hamiltonian, appropriate to the finite R4- 
space, is that described in Section 2.4. However, if one applies the unitary transformation to 
one obtains an effective Hamiltonian that is an A-body operator. To avoid this complication, an 
approximation is made in the NCSM. It consists in first finding only a two-body effective interaction 
v}j' eS ^ which is then used to replace the interaction term Vij of Eq. (94). The effective interaction V^ 2 ' eff ' 

is obtained by applying the decoupling condition of Eq. (84) to a two-nucleon Hamiltonian that 
arises from by restricting the sums to two nucleons only, keeping however the original mass number 
A in the interaction term: 
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The effective Hamiltonian is determined in a subspace of P, the P 2 -space (P 2 + Q2 = h), via the 
two-body transformation operator = <5 2 5 , ' 2 'P 2 . Then the two-body effective interaction is obtained 
by subtracting from the two-body HO Hamiltonian, i.e. 

T/[2,eff] _ tt[2] tt[2] 



(97) 



Figure 2: The various P and Q spaces relevant for the construction of the two-body effective interaction 
(see text). 



The corresponding V^ 2 ' eff ' is then used in Eq. (94). As is clear from Ref. [99], this procedure is equivalent 
to the case where (i) the similarity operator S of Section 2.4 is limited to a two-body operator 
and (ii) the effective Hamiltonian from the BCH expansion in Eq. (82) is truncated at the two-body 
operator level. This is different from the CCSD approach, described in Section 2.4.3. There one has 
S = S'M+jS'' 2 ', without a further truncation of the effective Hamiltonian. Due to the effective interaction 
the NCSM is no longer a variational method. The n > 2 terms neglected in S, as well as the n- 
body terms neglected in the corresponding effective Hamiltonian, could either increase or decrease the 
binding energy. On the other hand, as the basis space is increased the result converges to the exact 
solution. This becomes clear from the illustration in Fig. 2. At each Pa, the unitary transformation 
transfers information from a part of the Qa space, namely the Q2 space, to the P2 space. Thus E$ 
converges much faster when performing a calculation in P^-space than in a corresponding calculation 
with the bare interaction. For a sufficiently large Pa, covering practically the whole Hilbert space, the 
effective interaction coincides with the bare one, and one has the exact result. ^From the figure it is 
also clear that the same exact result could be achieved by enlarging horizontally the space where the 
unitary transformation operates, i.e. considering terms (as well as n-body terms in the effective 
Hamiltonian) with increasing n up to n = A. However, this is much more difficult, since one would 
need to know the entire n-body spectrum to construct the n-body effective interaction. It is evident 
that if three-body forces are present in the original Hamiltonian it is expedient to apply the procedure 
at least up to n = 3 [107]. 

Here we would like to mention another version of the NCSM, where the A-body basis is the trans- 
lationally invariant HO basis. Its application, however, is presently restricted to A = 3, 4. This version 
starts from the observation that in Eq. (94) the CM dependence is contained, as a separate HO Hamil- 
tonian, only in the one-body HO term. Therefore, the use of a translationally invariant basis expressed 
in terms of (A — 1) Jacobi coordinates automatically removes the CM contribution to the energy. The 
natural basis to use is the HO basis depending on the (A — 1) Jacobi coordinates, Ylf -1 
supplemented by the spin and isospin parts and totally antisymmetrized. 

The procedure of Suzuki [101] discussed in the previous section, which allows one to obtain the 
Hermitian two-body effective Hamiltonian, is applied to 
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on 



where ir± is the conjugate momentum of rj 1 . By using the SL transformation for one generates a 
two-body effective interaction 

V^ = m-(^\mnW) (99) 

that replaces Vjk in Eq. (94). 

Alternatively, and distinct from the usual NCSM approach, one could also obtain an effective po- 
tential V^*^ without an explicit Q dependence as follows 

vlr m = m _^i + ^JA^l nW ^. (100) 

This Jacobi coordinate formulation, implemented in coupled J-scheme, is more advantageous for A < 4 
nuclei. The single-particle formulation described above and implemented in M-scheme, is more ad- 
vantageous for A > 4. In fact the Jacobi coordinates prevent the use of Slater determinants, and 
antisymmetrization becomes much more difficult. 

Fig. 3 illustrates the convergence pattern of the ground-state energy of 4 He, taken from Ref. [110]. 
The results are obtained with the translationally invariant version of the method and with a realistic 
potential. It is interesting to note not only the Q dependence of the convergence but also how rapid 
the convergence can be for some fl values. Of course all the results have to converge to the same value. 

Concluding this section we would like to mention that momentum space NN potentials can easily 
be treated by using a HO basis in momentum space (see e.g. [111]). A review on recent developments 
of the method can be found in [112]. 
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Figure 3: N max dependence of the 4 He ground- 
state energy for a NCSM calculation with the CD- 
Bonn NN potential [109] and with various values 
of htt (from Ref. [110]). 
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Figure 4: K max dependence of the 4 He ground- 
state energy for a EIHH calculation with the MTI- 
III [182] (a) and AV14 [38] (b) NN potentials 
(from Ref. [114]). 



2.4.2 The Effective Interaction for Hyperspherical Harmonics (EIHH) 

The EIHH method, inspired by the NCSM approach, was introduced by Barnea and the authors of the 
present review article in Refs. [113, 114]. In this method one uses the HH basis to construct an effective 
interaction by applying the SL transformation to the Hamiltonian (see Section 2.4). 



on 



The idea is very similar to that of the NCSM approach, but with the following differences: (i) the 
HH basis instead of the HO basis is used, (ii) the P-space is defined by a maximal value K mBX of the 
grand angular quantum number (N = A — 1), and (iii) the EIHH two-body Hamiltonian that 

f2l 

undergoes the similarity transformation, H A , is a quasi two-body Hamiltonian, because it contains 
information about the dynamics of the entire A-body system. This is not the case for the HO two-body 
Hamiltonian, where the dependence on A in Vy of Eq. (95) is merely fictitious, since it does not add 
any dynamical information about the A-body system. 
The HH quasi two-body Hamiltonian is given by 

H [ 2 ] (p) = T K {p) + V A ,A-i(V2psm(a N )fj N ) , (101) 

where we make explicit only the spatial dependence of Va,a-\- The Jacobi coordinates are defined 
in normalized reversed order with r\ N = (ta-i — r^)/v^2, and T^-(p) is the collective hyperspherical 
kinetic energy of the entire A-body system (see Eqs. (55)). 

The unitary transformation is applied separately for each value of the hyperradius p, hence the ef- 
fective interaction becomes a function of p. In addition to its p-dependence, the HH effective interaction 
also depends on some quantum numbers of the residual system. For the more complicated case of a 
noncentral force we refer to Ref. [114] for an explicit expression, while here we illustrate the simpler 
case of a central potential model, where the matrix element of H A is given by 

([KnM\ P )\[K' n ]) = 5 [Kn][k , n] ^ Kn{Kn + 3N ~ 2) + , (102) 

(the various quantum numbers are defined in Section 2.3.1). One sees that the interaction depends 
on Kn-i and on L^-i , which are quantum number of the residual system, and on the quantum 
numbers Kn, Ln, and In, which characterize the entire A-body system. The various, just mentioned, 
additional many-body information contained in the HH effective interaction (obtained by subtracting 
the hypercentrifugal term) leads to a fast convergence of the HH expansion as can be seen in Fig. 4. 

Further improvements of the HH effective interaction have been made in Ref. [115] by introducing 
(i) a three-body effective interaction, and (ii) by including the hyperradial kinetic energy T p of Eq. (55) 
in the two-body Hamiltonian of Eq. (101). An extension of the HH effective interaction to non-local NN 
potentials has been carried out in Ref. [116], while an additional consideration of a 3NF is described in 
Ref. [117]. 

2.4.3 The Coupled Cluster (CC) Method 

The coupled cluster (CC) method originates from the field of nuclear physics [118] in the late 1950s. 
The approach was further developed during the 1970s. A summary of the status of the field in the late 
1970s can be found in Ref. [119]. Because of the difficulties posed by the strong short-range repulsion 
and by the tensor term in the nuclear interaction there were few nuclear physics applications then, see 
Ref. [120]. The most advanced of these were those of Mihaila and Heisenberg [121] who, at the end of 
the 1990s, were able to calculate the electron scattering elastic form factor of 16 O with realistic two- 
and three-body potentials. More successful, however, has been the use of the CC method in quantum 
chemistry [122], a field it was exported to 50 years ago and where it has had an enormous success. 
Today it has become a standard approach in quantum chemistry (for a review see [123]). Recently, 
starting with the work of Ref. [124], the CC theory has been reimported to nuclear physics. 

A first important observation is that in its original formulation the CC theory is not translationally 
invariant. Later, translationally invariant versions were formulated [120]. For clarity of presentation 
we will, in the following, first sketch the original idea and later comment on the problem of the CM 
spurious effects. 



One can understand the principal point of the CC method starting from what in quantum chemistry 
is called configuration interaction (CI) method. This is a variational expansion method. The variational 
wave function is expanded on the complete basis of z/-particle-z/-hole Slater determinants formed by 
Hartree-Fock or simply HO single-particle wave functions. In the language of second quantization the 
CI variational wave function is written as 

\y) N = (1 + C [1] + C [2] + C [3] + • • • + C [A] )\%) , (103) 

where |$o) is for example an A-body Slater determinant of the A lowest single-particle states. The 
i/-particle-z/-hole creation operator with v variables ijk... is defined by 

CM|$ >= ]T C%;;;v$t;;) (104) 

i>j>k..,a>b>c. 

with 

I^S:::) = alala\...a k a jai |$ ) , (105) 

where i < A and A < a < N. The restriction a < N means that the highest state that can be occupied 
is the N-th single-particle state. In principle only when iV — > oo is the full Hilbert space is covered. In 
practical applications one has only few terms in Eq. (103) and in most cases v is limited to two or 
three. Convergence both in v and N is very slow. 

In order to speed up the convergence for the ground state l^o) °ne uses the similarity transformation, 
of Section 2.4, applied to the Hamiltonian operator, 

E = (^o\H\^ ) = (^ \e s e~ s He s e- s \^o) = ($o|^|$o> , (106) 

where "H = e~ s H e s is the effective Hamiltonian whose mean value in the model space P, represented 
by !$())) coincides with E . As already discussed in Section 2.4, the most general form of the operator 
S is a sum of a one-body, a two-body up to an A-body operator, i.e. S = + + + • • • + S^ A \ 
In the language of second quantization the action of the various operators and on |$ ) is given 
by 

s [1] \$o) = E^I^) = E^>kl $ o) (107) 

i, a i,a 

S [2] \%) = E = E s^atala^o) , (108) 

i>j,a>b i>j,a>b 

and for with v > 2 accordingly. While in the CI the unknowns are the coefficients C^-, here they 
are the coefficients s^r;. There is clearly a relation among them. In fact, expanding the similarity 
transformation as (1 + S + S 2 /2\ + S 3 /3\ + ...) one can show that 

C [i] = 5 [i] > C® = S® + hsW) 2 , C^ = S^ + S^S^ + hs^) 3 (109) 

cm = s^ + ks [2] ) 2 + s^s^ + ks [1] )\s^) + -(s^)\ (110) 

where these equations can easily be extended to higher C^'. One notes that already a limitation of 
the expansion of S up to (CCSD) or to (CCSDT) induces correlations into \^o), which range 
from two-body (two-particle-two-hole) up to A-body correlations (A-particle-A-hole). This property 
is called size extensivity. However, from the equations above one also sees that a consideration of all 
operators up to order /i takes into account all correlations induced by operators with v < /i, 
but includes only a part of the correlations induced by operators with v > /i, namely those limited 



to the so-called disconnected components (the name refers to the graphic representation of the various 
correlations). 

In the CC method the problem is then to find the values of s"j^;;;. This is achieved by applying the 
decoupling condition QHP = 0, i.e by projection onto a sufficient number of excitations 

(^\U\%) = 0, ($g>|^|$ > = 0, ($g*|ft|$ o ) = (111) 

and further equations accordingly. Presently, in nuclear physics the CC calculations are limited to 
CCSDT, where the above three equations have to be considered. 

Here we would like to point out that it is expedient to use the BCH expansion of Eq. (82) in 
representing "H, since the expansion terminates after fourfold commutators when H has no more than 
two-body operators (after sixfold commutators when H has no more that three-body operators etc.). 

In order to illustrate the convergence rate and the role of doublets and triplets in a typical CC 
calculation, we show in Fig. 5 the results for the 4 He ground-state energy by Hagen et al. [125]. A 
few comments are in order here. Convergence is difficult to obtain for hard-core potentials, therefore 
in this example a very soft low-momentum interaction V\ owk [10] is used (the latter is obtained by 
implementing in momentum space the renormalization group evolution described in Section 2.4). Since 
a CCSDT calculation is at present prohibitively expensive compared to CCSD, a fifth order perturbative 
treatment of the triplets CCSD(T) has been implemented. It has been shown that for a specific choice 
of a soft potential, the addition of this contribution results in agreement with the FY precise ab initio 
result. Recently, a more sophisticated way of including the triplets, known as the A-CCSD(T), has been 
developed [127]. Performing such a A-CCSD(T) calculation [95] with the Idaho-N3LO (I-N3LO) [128] 
the 4 He binding energy came with a difference of about 0.15 MeV very close to the corresponding FY 
result. 

Now we turn to the problem of spurious CM effects in CC calculations. This is an old problem 
that was discussed already in the work of Zabolitzky [129], who used the common approach, already 
mentioned in Section 2.4.1, of adding a CM HO Hamiltonian to the intrinsic Hamiltonian. A differ- 
ent approach was used by Bishop et al. who formulated the problem and its solution in coordinate 
space [120]. Recently, it has been found numerically by Hagen, Papenbrock, and Dean [130] that in a 
sufficiently large model space, even when the model space is not a complete Nhfl space, the center-of 
mass wave function is approximately a Gaussian whose width varies little with the frequency Q of the 
underlying HO basis. The reported results on the 16 O nucleus open the door for a verifiable description 
of translationally invariant states for a large variety of model spaces and many-body methods. 

With the CC method it is possible to study excited states by making use of the equation-of-motion 
(EOM) formulation [131]. Essentially this consists in writing an excited state \^ k ) in terms of an 
excitation operator Q k acting on the ground state, i.e. \^ k ) = ^fcl^o)- By starting from the following 
two Schrodinger equations, 

H\* ) = £ |*o}, ( 112 ) 
H\^ k ) = £A|*o>, (113) 

and subtracting the former from the latter one obtains the EOM expression 

[H,n k ]\* ) = (E k -E )Sl k \V ). (114) 

In applying the EOM idea to CC theory, one starts from the CC ground state and writes the operator 
il k as a linear combination of the z/-particle-z/-hole operators (see Eqs.( 107), (108)). Therefore since 
Q k commutes with S one is then allowed to write 

[H,n k ]\$ ) = (E k -E )n k \$ ), (us) 
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Figure 5: N max dependence of the 4 He ground- 
state energy for a CC calculation with a V\ ow k 
NN potential; shown are CCSD (dotted) and 
CCSD(T) (dashed) results, also given FY re- 
sult [126] (full) (from Ref. [125]). 
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Figure 6: Results for e(r)-e(r = 0) from a GFMC 
ground-state calculation for 4 He (squares), 6 He 
(circles) and 7 Li (diamonds) with the AV18 NN 
potential [96] and the Urbana IX 3NF [139] (from 
Ref. [139]). 



or, equivalently, 

nn k \%) = E k n k \$ ) . (116) 

Therefore one can construct the H matrix on the z/-particle-z/-hole basis and then diagonalize it to obtain 
the excitation spectrum. Special care has to be paid in the last two steps, since T-L is non- Hermit ian. 



2.5 A Formulation Suitable for Monte Carlo (MC) Approaches 

In this section we describe four methods for which the formulation of the quantum mechanical many- 
body problem lends itself to a stochastic approach. These methods, which are generally referred to 
as Monte Carlo techniques, are the Green's function Monte Carlo (GFMC), the chiral effective field 
theory on a lattice (LCEFT), the Monte Carlo shell model diagonalization (MCSMD), and the varia- 
tional Monte Carlo (VMC). The GFMC method is based on the path integral formulation of quantum 
mechanics, as is the diffusion Monte Carlo (DMC). The differences between the GFMC and DMC 
methods are actually very small so that in the literature they are often interchanged. The LCEFT is 
essentially a DMC approach as well, except that the dynamical degrees of freedom are fields rather 
than particles. The GFMC and LCEFT methods are based on the Euclidean time (imaginary time) 
evolution of the system. On the contrary the MCSMD, although inspired by the imaginary time formu- 
lation, is still effectively a variational method. In fact MCSMD starts formally from an imaginary time 
evolved trial function, which after some manipulation, leads to an expression that suggests a way of 
constructing a variational SM basis. The imaginary time becomes just one of the non-linear parameters. 
The fourth approach, VMC, is a fully variational method, where the MC technique is used to evaluate 
the many-dimensional energy functional integrals. The variational wave function that is obtained with 
this method is very important since it serves as starting trial function for the GFMC imaginary time 
evolution. 



2.5.1 The Green's Function Monte Carlo (GFMC) Method 

The first application of a GFMC idea is due to Kalos [132] who in 1962 calculated the ground state 
of three- and four-body nuclei using some simple central interactions like square well and Gaussian 
potentials. The application of the GFMC algorithm to nuclear Hamiltonians containing spin and 
isospin operators was proposed in 1987 by Carlson [133]. To deal with a large number of nucleons and 
Hamiltonians containing spin and isospin operators the auxiliary field diffusion Monte Carlo (AFDMC) 
method, which uses the Hubbard- Stratonovich (HS) transformation [134, 135] to linearize quadratic 
terms in the evolution operator, was introduced by Schmidt and Fantoni [136]. 

With the GFMC method it is possible to aim directly at the calculation of ground-state and low-lying 
energies. In fact it is shown below that the ground-state energy E can be obtained by 

E = lim e(r) , (117) 

T— S-OO 

where r is the so-called "imaginary time" defined as r = it and 

e(T)= <¥|¥(r)> • (118) 

Here |\&(t)) represents the evolution in imaginary time of a function \^/) (in principle any function, 
called a trial function) 

|^r(r)) = Ne- iHt \V) = Ne~ HT \^) , (119) 

where the normalization N is chosen as N = e rE °. Therefore, instead of solving the Schrodinger equation 
to calculate the ground-state energy of a system, one determines E via 

( s S>\He~( H ~ E °"> T \^\ 

The proof for this assertion is rather simple. For this purpose one formally writes the trial function 
as a linear combination of the eigenfunctions |^ 7 ) of H: 

|*> = |>c 7 |* 7 >. (121) 

Inserting this expression in Eq. (118) one finds 



2 o -(£ 7 -£ )t 



e(r) = ^ = ^-j . (122) 



When the imaginary time tends to infinity and under the condition that the trial function is not 
orthogonal to the ground state, the only contribution to the sum/integral that survives is the ground- 
state energy: 

e(r^oo) = £ - (123) 

The reason why this formulation takes the name of Green's Function MC is clear if one represents in 
configuration space the function |^(r)), necessary to calculate e(r) numerically. In fact in configuration 
space one has 

|tf(r)) = e ^ H ~ Eo)r \^) 

= J J dRdR'\R)(R\e- {H - Eo)T \R')(R'\y) = J J dRdR'\R)G{R, R', r)(R'|*) , (124) 



or. 



where R represents all the particle coordinates and G(R, R', r) is just the Green's function of the 
imaginary time evolution operator. 

In order to evaluate |^(r)) one writes the evolution operator as 



e -(H-E )r 



e -(H-E )Ar 



(125) 



and therefore the Green's function becomes 



G(R,R,r) = J ••• J c£R n c£R n _i • • • dR 1 G f (R, R n , Ar)G(R n , R n _i, At) ■ ■ ■ G(K 2 , Ri, At)G(Ri, R', At) . 

(126) 

At each time step the evolution is implemented by the propagation of a finite collection of points 
((^-functions) approximating the function in R. 

Dividing the evolution into small imaginary time steps is crucial in order to write the explicit 
expression for G(R, R n , At) as [137]. 



1 \¥ f mfR-R^ 2 



G(R, R', At) = — — exp \^ 1 - - (V(R) + V(R')) - E At . (127) 
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Such an expression is correct only to the order At. In order to better understand the above expression 
one can consider the Schrodinger equation in imaginary time as equivalent to a diffusion equation with 
an absorption term governed by the Hamiltonian. In the zero mass limit, i.e. when the kinetic energy 
becomes dominant, the solution of this diffusion equation is 

*(R, t)= j dR'G(R - R', r)*(R', 0) (128) 

with the Green's function given by 
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On the other hand, the consideration of the infinite mass limit gives the term [(V(R) + V(R'))/2 — Eq] 
proportional to At in the argument of the Gaussian. It is clear that the sum in the exponent of a 
kinetic and a potential term is exact only for small increments of the imaginary time, since the two 
operators do not generally commute. 

There are a few delicate points in the numerical application of the GFMC procedure. One of them 
has to do with the choice of the normalization constant iV = e rE °, which is just a function of the 
ground-state energy itself. Such a choice is crucial to keep the denominator in Eq. (118) equal to a 
finite value (actually equal to one) in the infinite r limit. So one of the delicate points consists in 
keeping the normalization under control without knowing the value of Eq. 

Another important point is the so-called "sign-problem" that arises in the case of fermions. In fact 
the time evolution leads in principle to the state of minimal energy, which is the bosonic ground state. 
There are various ways to cure this problem approximately. The most widespread is the fixed node 
approximation. It is based on the fact that an antisymmetric continuous wave function must become 
zero on some hypersurface of dimension (3N-1). Since one does not know this hypersurface, i.e. one does 
not know the correct nodal structure of the exact fermion wave function, one imposes that the evolved 
function has the same nodal structure as the antisymmetric trial function. The error made by this 
approximation is not known. One can only prove that the value of the energy obtained starting from 
a variational trial function (like e.g. the VMC discussed in Section 2.5.4) lies between the variational 
result and the exact one. We refer to Ref. [138] for a comprehensive discussion of the sign problem. 

The GFMC method requires a very good trial function as a starting point, in order to be able to 
reach convergence, while increasing the imaginary time. Therefore the VMC result is often used as a 
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trial function. To have an idea of the quality of the results for E in the r infinity limit and of the rate 
of convergence in Fig. 6 we show one of the results from Ref. [139]. As expected the errors grow with r 
and A. It is interesting to note that the 4 He points do not show visible errors. This is due to the fact 
that the sign problem is avoided in this case. In fact for 4 He the antisymmetry of the wave function 
resides entirely in the spin-isospin Slater determinant. Therefore the spatial part is symmetric and the 
problem is equivalent to a boson problem. 

The GFMC technique can be applied to excited states as well. If they have quantum numbers 
different from the ground state, the procedure is exactly as described above, provided that the trial 
function has the desired excited state quantum numbers. It is also possible to treat a few excited states 
with the same quantum numbers as the ground state or any other reference state [140]. In this case one 
proceeds as follows. One first constructs a r-dependent Hamiltonian matrix H(t) between normalized 
states 

' ' l*.(V2)||*,(r/2)| 1 

This is possible because in GFMC one can compute mixed expectation values starting from states 
and |\F,) obtained by diagonalizing a VMC Hamiltonian and propagated separately. One then solves the 
generalized eigenvalue problem for H(t). For large values of r one recovers the exact discrete spectrum 
and a discretized continuum spectrum. 

A different technique is largely used in atomic, molecular and cluster physics. Here one starts from 
the imaginary time evolution of a state which is prepared applying a projection operator n on an initial 
trial function that selects the proper quantum numbers. More explicitly, what is calculated by GFMC 
for many values of r is the quantity 

£(r) = (^|nV^-^n|^). (i3i) 



It can be easily shown that C(t) is equivalent to the Laplace transform of the response function R(E) 
of a system, which is initially in a state ipi, to a perturbation induced by n (see Eqs. (140) and (144)). 
This means that knowing R(E) one also knows the excitation spectrum of H induced by n. This is 
recovered from its Laplace transform by means of different algorithms (see e.g. [141]). Due to the ill- 
posed nature of the inverse Laplace transform, however, these algorithms may have a scarce efficiency 
(see also Section 2.6.1). 

The GFMC technique can also be applied to low-energy scattering, where only a single two-fragment 
channel is open. As described in Ref. [142] one proceeds as follows. First a VMC calculation is carried 
out for a given scattering energy E in such a way that at a large distance, Rq, the wave function has 
the typical asymptotic form consisting of internal wave functions for the two fragments and a relative 
wave function for a given angular momentum, where the radial part is a proper superposition of Bessel 
and Neumann functions. This leads to a well-defined logarithmic derivative 70 at Rq. In a subsequent 
GFMC calculation the VMC wave function is used as trial function. The GFMC simulation is restricted 
to a volume corresponding to cluster separations less than R and the logarithmic derivative 70 is used 
as boundary condition. The GFMC procedure then leads to the proper energy E corresponding to the 
boundary condition. Having determined Rq, 7o, and E the phase shift is uniquely defined. 

The applicability of the method is limited to local potentials in configuration space. Potentials 
in momentum space, and in general non-local potentials, can generate great problems in finding an 
expression for the propagator G(FL, R n , At) that can be sampled. Therefore phenomenological local 
configuration space potentials have been constructed ad hoc for application of the GFMC approach to 
nuclear systems; the modern realistic AVI 8 potential is one of them. 
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2.5.2 Lattice Chiral Effective Field Theory (LCEFT) 



Effective field theory provides a systematic approach to studying low-energy phenomena in few- (and 
many-body) systems. Like in all effective field theories, there is a wide disparity between the long- 
distance scale of low-energy phenomena and the short- distance scale of the underlying interaction. One 
can characterize the low-energy phenomenology according to exact and approximate symmetries and 
low-order interactions according to some hierarchy of power counting. In the case of nuclear physics 
and for nucleon momenta comparable to the pion mass, both the contribution from nucleons and pions 
must be included in the effective theory. 

Recently an innovative method that combines the MC formulation with the chiral effective field 
theory of the nuclear force has been introduced (for a review see Ref. [143]). The approach is similar to 
lattice QCD, but the fundamental fields are replaced by effective nucleon and pion fields. The nuclear 
interaction is described by an effective chiral potential and the inverse of the lattice spacing serves as 
an ultraviolet regulator of the theory. The ground-state energy of the nucleus is also here filtered out 
at large imaginary time. 

In practical applications one builds up the action S in terms of nucleon and pion fields. If the action 
is given by the functional S of field configurations, then the time ordered vacuum expectation value of 
a functional H is given by 

{H) = jv^m ' (132) 

where the symbol / D<\) represents the infinite-dimensional integral over all possible field configurations 
on the entire space-time. 

In view of an evaluation of (H) on the lattice one can express the time integral of the Lagrangian L 
in the action S as a sum of N small time intervals At, and therefore %S as a sum of N small Euclidean 
(imaginary) time intervals At: 

iS = lCAt k = E CAr k • (133) 

k k 

In a four-dimensional lattice evaluation of (H) the time interval Ar corresponds to the imaginary time 
lattice spacing a T . 

/,From the technical point of view one then proceeds in essentially the same way as in the GFMC 
(compare Eq. (132) with Eq. (120)). Due to the appearance of quadratic forms in the fields in the 
so-called transfer matrix, i.e. of the normal-ordered exponential of the Hamiltonian over one temporal 
lattice spacing, one makes use of auxiliary fields, like in the AFDMC. The main difference between the 
LCEFT approach and the GFMC (or AFDMC) consists in the different choice of the dynamical degrees 
of freedom: particles in the case of GFMC or AFDMC, fields in LCEFT. 

In principle potentials derived in the EFT framework could also be used in the GFMC method 
provided that their representation in configuration space is given. 



2.5.3 Monte Carlo Shell Model Diagonalization (MCSMD) Method 

There are two different methods that can be confused because they are both referred to as Monte Carlo 
shell model. The first, described in Ref. [144] is an extension of Eq. (117) to large SM problems with 
core. The second, MCSMD, introduced in Ref. [145] and illustrated in [146], is in some way similar to 
the SVM discussed in Section 2.3.3, in that it is based on a stochastic generation of the basis. 

The MCSMD method has been inspired by the AFDMC, but it does not make use of Eq. (117). 
Instead the method consists in creating a basis starting formally from an imaginary time evolved trial 
function (taken as a HO slater determinant |^o)) 

\$) = e-* H \* T ), (134) 
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and using the HS transformation [134, 135]. This linearizes the two-body interaction by auxiliary fields 
and allows to be rewritten as 

\ij) = J daN(a)e-™ h \V ) , (135) 

where the operator h is now a one-body operator. At this point, however, one does not proceed to the 
t — > oo limit, but instead one observes that the vectors 

\tf>(a)) = N(a)e-™ h \* ) , (136) 

could be good candidate states for a basis. In order to explain the MCSMD procedure, let us assume to 
have already chosen a certain number of basis states by random a values. The ground-state energy is 
calculated. Another vector is added to the basis by choosing a new random value of a, and after having 
orthogonalized it with respect to the other basis vectors the Hamiltonian is diagonalized again. Only if 
the ground-state energy shows a significant decrease the state is kept. The procedure is repeated up to 
convergence. During the MCSM generation of the basis vectors, symmetries, e.g. rotational and parity 
symmetry, are restored before the diagonalization. 

2.5.4 The Variational Monte Carlo (VMC) Method 

Another variational method that uses the stochastic approach is the variational Monte Carlo (see 
e.g. Ref. [147]). In this case the trial function, written in configuration space, consists of a sum of 
differently correlated antisymmetric single-particle states. Correlations of a various nature are obtained 
by operating with a certain number of two- and three-body correlation operators. They reflect the 
operator form of the potential and have an influence on the short-range part of the wave function. 
Appropriate boundary conditions are imposed at long range. The stochastic nature of the approach lies 
in the Metropolis Monte Carlo evaluation of the energy functional multi-dimensional integral. 
More specifically the variational trial function has the form 

*j> • (137) 

The Uij and are non-commuting two- and three-nucleon correlation operators 

^,- = E«p(^)Q?-, (138) 
v 

where the Ofj are the same operators that appear in the NN potential used. The Wijk have the 
spin-isospin structure of the dominant parts of the 3NF, as suggested by perturbation theory, S is a 
symmetrization operator and ^ j is the fully antisymmetric wave function that has already been Jastrow 
correlated. For example, for s-shell nuclei one writes 

<5> A (JMTT 3 )) , (139) 

where \&a(JMTT 3 )) is a simple spin-isospin Slater determinant and the various / are correlation 
functions. 

For a generic nucleus \^fj) is more complicated because it is a linear combination of N channel 
basis functions characterized by total orbital angular momentum L, total spin S and principal quantum 
number set [n\. 

The overall wave function is translationally invariant. This is clear when the single-particle basis 
is a simple spin-isospin Slater determinant, or a linear combination of them, since the correlation 
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Table 7: 4 He ground-state expectation values in units of MeV for ground-state energy E, kinetic energy 
T, Vij (AV14 NN potential [38]), and V ijk (Urbana VIII 3NF [148]) from VMC [149] and GFMC [90] 
calculations (GFMC result as given in [149]). 





E 


T 




Vijk 


VMC 


-27.2(2) 


106.6(8) 


-129.7(7) 


-4.84(9) 


GFMC 


-28.3(2) 


113.3(20) 


-136.5(20) 


-5.8(3) 



functions are functions of distances between particles. However, even when the single-particle basis 
is a Slater determinant of single-particle wave functions from a parametrized one-body potential, one 
can implement translational invariance numerically. In fact at each step of the Metropolis sampling 
the values of the A particle positions are substituted by the positions with respect to the CM 
— Vi — 1/A J2i r i an d the kinetic energy is renormalized in an analogous way. 

In Table 7 we show VMC results for various 4 He ground-state quantities in comparison to GFMC 
results. One notes that the ground-state energy underestimates the precise ab initio result by only 1 
MeV, whereas the differences amount to about ±7 MeV for expectation values of kinetic energy and 
NN interaction. The results for the 3NF differ by 1 MeV. Summarizing one can conclude that the VMC 
is already rather close to the precise ab initio results and thus supplies a very good trial wave function 
for the GFMC calculation. 

2.6 Bound-State Formulations of the Continuum Problem 

Various formalisms discussed in the previous sections have been applied in the past to ab initio cal- 
culations of nuclear observables for nuclei with A > 3. However, only energies below the threshold 
for three-body break-up were considered. The solution of the quantum mechanical problem in the far 
continuum i.e. where more than two fragments exist is much more complicated, since one then has to 
deal with the full many-body scattering problem. Several methods have been proposed to reduce the 
continuum state problem to a bound state one [150, 151, 6, 152], which have in common the use of a 
complex expression for the resolvent operator. In Refs. [152, 153] a numerical extrapolation to the real 
axis is performed. In the following we will discuss the two methods [6, 154, 155], where one remains on 
the complex plane. They enable calculations of perturbation induced reactions to the continuum and 
at the same time avoid an explicit determination of continuum wave functions. 

2.6.1 The Lorentz Integral Transform (LIT) Method 

The LIT method was first formulated for perturbation induced reactions involving transitions to the 
continuum [6]. The success of the method has been such that it has lead to ab initio results for 
observables of A = 4, 6, 7 nucleon systems beyond the three-body break-up threshold and even in the 
full continuum. 

To explain the LIT method consider a typical observable for a perturbatively-induced reaction, 
namely the inclusive response function 

R(E) = \(yB(E y )\OMEo))\ 2 6(E y - E ) , (140) 

where O is an operator inducing a transition from the ground state *&(E ) with ground-state energy 
Eq to a state ^(Ej) with energy E 1 . Performing an integral transform with a well-defined and smooth 
kernel K(a, E) one has 

$(cr) = J dE R(E) K(a, E) (141) 

a n 



which by inserting the relation (140) becomes 



$(a) = ^d 1 {^{E Q )\0^{E 1 ))K{a,E)(^{E 1 )\0\^{E G )) 

= ^^{E^K^H^iE^^iE^O^iEv)). (142) 

Use of the closure property of the eigenstates of H, 

£di\*(E 1 )){*(E 1 )\ = l, (143) 

then leads to 

$(<t) = (*(E )\K(<t,H)\*(Eo)). (144) 

The principal idea of the integral transform method consists in the calculation of the transform $(<r) 
without knowing R(E) and next in the inversion of the transform in order to obtain it. For this purpose 
one needs an appropriate kernel K(a, E). In the past the Stieltjes transform with K(a, E) = 1/(E + a) 
with a > 0, proposed by Efros [151], and the Laplace transform with K(a,E) = exp(—(rE), used in 
Ref. [156], were considered in the past for nuclear physics problems first. Both transforms suffer a 
drawback in that information about the response R at a given energy E is dispersed over a rather large 
a range. This makes the inversion of the transform rather unstable as discussed in detail in Ref. [157] 
for the Stieltjes transform. 

Considerable improvement of the inversion results from kernels that distribute the strength over 
only a relatively small a range. Such a request is satisfied by kernels that are representations of the 
5-function. Therefore Efros and the authors of the present review article have proposed the Lorentz 
kernel [6] 
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K ( a > ^ = {E - E - a R f + a] = E - E - a R - taj E-E -a R + taj ' (145) 

where 07 is a constant chosen to be sufficiently small and a R > 0. Insertion of this kernel into Eq. (144) 
leads to the following result for the LIT: 

L(a R ,a I ) = (*\V). (146) 

The localized LIT function ^ is obtained from the LIT equation: 

(H-E -<T R -i<T I )\V)=O\*(E )). (147) 

The LIT equation represents a bound-state like equation and thus can be solved by applying one of 
the usual methods for a solution of a bound-state problem. Once L(a R , 07) is calculated R(E) can be 
obtained by inverting Eq. (141). The reduction from an originally continuum-state to a bound-state 
like problem is an enormous simplification and allows the calculation of electroweak observables which 
otherwise, until the present, cannot be calculated. 

A benchmark test for the triton total photoabsorption cross section in unretarded dipole approxi- 
mation has been performed in Ref. [158]. For this purpose the transform L(a R , 07) has been calculated 
for a R — 0, 1, 2, 300 MeV and at a fixed 07 of 10 MeV. This allowed the integral transform to 
be inverted with great accuracy so that the unretarded dipole response function i? Dunrct (_E 7 ) could be 
determined with high precision. Finally by multiplying i?D iUnret (i? 7 ) with a trivial kinematical factor 
one obtains the total photoabsorption cross section in unretarded dipole approximation. In Fig. 7 the 
small inversion error range is represented by the space between the two almost overlapping curves, while 
the dots represent the Faddeev result calculated at various energies. A nice agreement between the LIT 
and Faddeev results can be observed. 
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Figure 7: Results for the total 3 H photoabsorption cross section with the AV18 potential from LIT/CHH 
(full) and Faddeev (dots) calculations in unretarded dipole approximation (the two curves represent the 
uncertainty of the inversion of the LIT); from Ref. [158]. 



The LIT approach can also be applied to exclusive processes, for which details we refer to Ref. [66] . 
Applications have been made for the reactions 4 He(7,N) 3 X [159], 4 He(e,e'p) 3 H [160], and 4 He(e, e'd) 2 H 
[161]. The applicability of the LIT to photon scattering has also been established in Ref. [162]. 

A remark concerning the inversion of the LIT is in order here. The inversion represents a so-called 
ill-posed problem. The mathematical term ill-posed is somewhat misleading, since ill-posed problems 
can be solved reliably. For a detailed discussion of this issue we refer to Ref. [163], where it is pointed 
out that the LIT method has to be understood as an approach with a controlled resolution. 

The solution of Eq. (147) can be obtained by matrix inversion after expanding the function |\&) on 
localized functions. However, in practical calculations one can avoid the matrix inversion. In fact one 
can rewrite L(<tr, 07) as 

L(a R ,aj) = —Im{($(Eo)\Oi 1 , O\*(E ))\ , (148) 

07 I PR + E — tl + lOj J 

and apply the Lanczos algorithm with pivot 0\^(Eq)) [164]. 

It is interesting to note that if the Lorentz kernel tends to the 5-function (07 — > 0, a R + E = E), 
Eq. (148) corresponds to the relation between R(E) and the so-called linear response x(E), 

R{E) = _I/ m{x (£) } = -l/ m ((^)|Ot 1 O\*(E ))\ . (149) 

71 71 i Ej — ti + ie ) 

We would like to stress that in the LIT method the continuum problem reduces to a bound-state like 
problem, due to the fact that the resolvent 1/(E — H + ie) acquires a finite imaginary part in the 
denominator. 

2.6.2 The Complex Scaling Method (CSM) 

Another approach whose goal is the rigorous reduction of continuum problems to bound-state-like 
problems is the complex scaling method. This method appears to have been first introduced in 1944 
by D.R. Hartree et al. for solving problems relevant to air defense research and development [165]. It 
was rediscovered in the late sixties/early seventies [150, 166, 167] to study scattering amplitudes and 
resonances and is still an active field of current research in atomic and molecular physics (for a review 
see Ref. [168]). 

An 




The method is based on a similarity transformation performed by the complex scaling operator 

S(0) = e ierd/dr (150) 



such that 

S(9)f(r) = f(re* d ) (151) 

for any analytical function f(r). This corresponds to a transformation of r and p operators such that 
r — > e l9 r and p — > e l9 p. 

The complex scaling is particularly suited for studying resonances since it associates the resonance 
phenomenon with the discrete part of the spectrum of the complex-scaled Hamiltonian. (The complex 
eigenvalues do not depend on 9, although the resonance complex-scaled eigenfunctions do). 

In nuclear physics ab initio calculations the CSM has only been applied in Ref. [169] for calculating 
elastic and three-particle break-up observables in the 3N system, and in Ref. [170] for investigating 
four-neutron resonances. Here, however, we would like to mention an interesting application of the 
complex scaling idea to the direct calculation of R(E). 

In Ref. [155] the complex scaling transformation has been used in the expression of R(E), 

R(E) = -I/ m |(*(^)|Ot^ (e) __l__ £(0)(9|*(£o)>} , (152) 

where H{9) = S(9)HS^(9). In this way the resolvent acquires a finite imaginary part. Therefore, within 
a suitable range of positive 9, one is allowed to solve the eigenvalue problem for H{9) by an expansion 
on square integrable functions. If sharp resonances exist, 9 has to cover the resonance poles on the 
complex energy plane [168]. 

A very nice feature of the CSM calculation of R(E) is that one has a direct access to it and no 
inversion is required. On the other hand, while in the LIT case the stability of the inversion has to be 
studied, in the CSM the stability of R(E) with respect to 9 has to be examined, 

The CSM approach appears to be very promising for calculating perturbative inclusive reactions 
and deserves to be investigated further. 



3 Discussion of Results 

In this section we present an overview of results obtained from ab initio calculations for nuclear systems 
with A > 4. A complete overview of the field is beyond our scope. Instead, here our aim is, on the one 
hand to illustrate the development of the field from a historical perspective and on the other hand to 
assess today's state of the art, paying particular attention to benchmark calculations. 

It goes without saying that modern-day methods and techniques in few-body calculations could 
not have been developed without the continuous progress made over, at least, the last fifty years. 
Therefore we think that, in addition to reporting more recent results, we have also to pay tribute to 
those researchers who have made the first attempts to solve the nuclear many-body problem with ab 
initio calculations using the most realistic potentials of the time, thus laying out the basis for further 
developments. 

Here we are dealing with nuclear systems with A > 4, but one should not forget that progress in 
treating the 3N system had an important influence on the development of the field of more complex 
few-body nuclei (for reviews see Refs. [23, 171, 172]). 

One should note that the spectra of nuclear systems with A > 4 are generally much richer than 
those of the 2N and 3N systems. Let us take as example the 4N system with isospin T = 0: One has the 
4 He bound state and five different continuum channels with five different thresholds, namely two (3+1) 
channels ( 3 H-p, 3 He-n), one (2+2) channel (d-d), one (2+1+1) channel (d-np), and a (1+1+1+1) 
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channel (2n2p). This can be compared to the 3N system in the isospin T = 1/2 channel with projection 
T z = 1/2, where one has one bound state ( 3 He), but only two different continuum channels ((dp) and 
(npp)). The comparison shows that the difference between the 3N and 4N systems does not just consist 
in one additional nucleon, but also in a more complex structure of the excitation spectrum. 

Our discussion will show that whereas enormous progress has been made over the years in few- 
nucleon calculations, still quite some work remains to be done. The most visible advancement has been 
obtained for bound-state calculations, where even A > 10 systems have been tackled. On the other 
hand continuum-state problems are considerably more complicated. For example, with the exception 
of the LIT method, common ab initio methods have not yet succeeded in controlling the complete 4N 
continuum up to pion threshold. Rather, the energy has been limited to the range below the three-body 
break-up threshold. Thus, it is not surprising that also for the A > 4 nuclear systems such ab initio 
methods are presently only able to treat the energetically lowest two-fragment break-up channels. 

The following discussion is organized as follows. In Sections 3.1 we describe the main aspects in the 
historical development of A = 4 calculations for the time period up to the late 1980s when a new phase 
in the 4N calculations had its beginnings. This modern period is then discussed by considering first 
bound-state calculations in Section 3.2.1 and then continuum-state calculations in Section 3.2.2, where 
we will distinguish between hadronic and perturbation-induced reactions. 

3.1 4N Problem: The Early Phase 

Probably the beginnings of serious a particle bound state calculations using semi-realistic NN potentials 
can be placed near the end of the 1960s. Variational calculations [173, 174, 175] were performed with 
local central NN potentials, with and without spin dependence (some of them also fit to describe the 
NN S- wave phase shifts). The trial functions were constructed using Jastrow type correlation functions. 

In the following years non-local separable potential models became more popular. In fact, for the 
two-body problem, it was long known [176] that non-local separable NN potential models, like the 
Yamaguchi NN potential models [176], reduce the problem to an algebraic equation. Somewhat later it 
was also realized (see e.g. Ref. [177]) that the three-body problem leads to an effective two-body problem 
by use of separable NN potentials. The generalization of such a procedure, namely the reduction of the 
problem from A to A - 1, was pointed out for the AGS formalism in Ref. [42] (see also the discussion in 
Section 2.2). The equations obtained for the (A-l)-system then led to quasi-potentials, which again 
can be described by an expansion in separable terms. As was shown in Ref. [42] a repeated application 
of such a procedure increasingly reduces the complexity of the equations. 

A first 4N calculation along these lines was performed by Alt, Grassberger, and Sandhas in Ref. [178] 
who used Yamaguchi potentials for the 3 Si and 1 s NN partial waves and took separable approximations 
for the (3+1) and (2+2) subsystem amplitudes. The resulting AGS equations were however solved only 
rather approximately. Similar reductions with separable models were also made for the FY equations 
by Kharchenko and Kuzmichev [179] and by Narodetsky, Galpern, and Lyakhovitsky [180]. Spin and 
isospin degrees of freedom were included in the work of Ref. [180] but not in that of Ref. [179]. The 
corresponding integral equations were solved with the Yamuguchi potential as a NN interaction, and by 
further assuming that the (3+1) and (2+2) subamplitudes are dominated by s-waves. In addition, the 
convergence of the calculated observables with respect to the separable expansion of the subamplitudes 
was checked. Shortly after, Tjon [181] performed essentially the same calculation as in Ref. [180], but 
used instead of a separable NN potential, a local s-wave potential of Yukawa type (MTI-III potential 
[182]) represented by a separable expansion. A rather realistic result of 29.6 MeV for the 4 He binding 
energy was obtained. In addition, by using various approximations to the local interaction, Tjon could 
show that there was a linear relationship between the binding energies of 3 H and 4 He. This was later 
confirmed for many other interactions and the linear relationship received the name Tjon-line. On the 
other hand it was not yet clear whether the separable ansatz for the amplitudes of the subsystems was 
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Table 8: 4 He binding energy in MeV for the MTV potential [182] with various methods as published 
in Ref. [199] in 1982. 



Method 


GFMC [198] 


VMC [200] 


CC [129]) 


ATMS [204, 205] 




31.3(2) 


31.19(5) 


31.24 


31.3 



reliable. This was investigated by Gibson and Lehmann [45] by directly solving the coupled two-variable 
integral equations resulting from the Schrodinger equation for the Yamaguchi NN potential, and, in fact 
the results of Ref. [180] could be reproduced. 

Here we would like to point out that more realistic NN potentials models than those mentioned 
above were also in use. Demin, Pokrovsky, and Efros [183] applied the HH method (potential harmonics 
approximation for K > 6) and obtained a rather good convergence of the expansion with the EH [184] 
and GPT [185] potentials, which led to binding energies of 23.03 and 26.76 MeV, respectively. This 
range of binding energies for realistic NN models is confirmed in GFMC and Faddeev calculations (see 
below and also Table 2). Applying the CC technique, Zabolitzky presented a calculation [186] with the 
Reid potential [187] with a full and partial consideration of the potential in the two-body (CCSD) and 
three-body (CCSDT) parts, respectively. A 4 He binding energy of about 21 MeV, and, in an improved 
calculation, of 23.7 MeV [119], were obtained. It should be mentioned, however, that Zabolitzky's 
approach was not completely translationally invariant for parts of the calculation. 

With time FY calculations became somewhat more realistic. In the late 1970s tensor force effects 
were included in an approximate way by Tjon [188]. He calculated the 4 He binding energy with a 
truncated Reid potential retaining in the NN t-matrix only the s-wave amplitudes, but such that the 
3 s\- 3 di tensor force was included perturbatively (so-called too approximation). This led to a binding 
energy of 19.5 MeV. Shortly after, the relative good quality of such a truncated t-matrix approximation 
was shown for the 3N case by Gibson and Lehman [46], and there, in addition, the same approximation 
was used to calculate 4 He binding energies with separable NN potential models. 

One focus in the early literature of FY/AGS 4N calculations concerned the best choice for the 
expansion method for the subsystem amplitudes. This was discussed in quite a few publications by 
various authors [189, 190, 191, 192, 193, 194]. Still another point of discussion concerned the necessity 
of introducing separable expansions not only for the (3+1), but also for the (2+2) subsystem. As men- 
tioned by Narodetsky, Galpern, and Lyakhovitsky [180] the expansion was made because of convenience 
and that an exact method, the convolution method, exists for the (2+2) case. An application of the 
convolution method is found in the field theoretic approach [195]. By applying the convolution method 
to the AGS formalism it was shown by Haberzettl and Sandhas [196] that the field theoretic formalism 
of [195] corresponds to the AGS formalism in the case that only single-term separable approximations 
for NN and (3+1) amplitudes are employed. 

In the early 1980s other groups appeared on the scene. HH calculations with central NN potential 
models were performed by Ballot [197]. At the same time the first more realistic GFMC calculations 
in nuclear physics appeared. These were carried out by Zabolitzky, Schmidt, and Kalos [198, 199], 
who considered 3N and 4N systems and employed the MTV potential [182] in addition to a simple 
central 3NF model [200]. Nearly concurrently variational calculations were initiated for 4 He from 
the Urbana-Argonne group. The trial functions were constructed variationally from various choices 
of pair correlation operators and/or of three-body correlation operators (see also Section 2.5.4). NN 
potentials, like MT models and the V8 version of the Reid potential, were considered without [201], 
and with the above mentioned simple 3NF [200]. The models for the NN potentials and for the 3NF 
were improved in subsequent variational calculations by Carlson, Pandharipande, and Wiringa [202] 
and by Schiavilla, Pandharipande, and Wiringa [203]. It is worthwhile mentioning that due to the 
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above mentioned developments a benchmark for the MTV potential was established in the early 1980s. 
The corresponding results are listed in Table 8. The various results agree quite well. Though the 
potential model is rather simple one may regard the agreement, obtained with four different methods, 
as a milestone in the history of the 4 He binding energy calculations. The advantage of such established 
precise ab initio results is self-evident: they serve as a standard for future calculations. Later the 
various methods that appear in Table 8 were further improved and applied to more realistic potential 
models. Only the variational ATMS method, to our knowledge, has never been applied to more realistic 
interactions. 

Improvements in AGS calculations were made by allowing more exact representations of the s-wave 
(3+1) amplitude [206]. Articles, in a common book, by Ferreira [207], Oryu [208], and Plessas [209] 
clarify the important role played by separable expansion techniques for both subamplitudes and NN 
potentials at that stage of the treatment of few-body problems. 

All the above mentioned FY/AGS calculations tackle the 4N problem in momentum space. Config- 
uration space FY calculations were initiated by Merkuriev, Yakovlev, and Gignoux [28]. They derived 
the differential equations for the FY components and, in addition, described their asymptotic form. 
Four-body binding energies were calculated with local s-wave NN potential models. A rather good 
agreement with results from other calculations [181, 197] was obtained. 

Important progress for the GFMC calculations was made by extending the technique to spin- 
dependent forces [133]. Calculations for potentials of V6 form (central potential with spin and isospin 
dependences plus tensor force) and simpler models were performed. A comparison with Faddeev cal- 
culations for the 3N system showed that precise results could be obtained with the GFMC technique, 
whereas variational results were not of such good quality. In addition the following results for the 4 He 
binding energy were found for the VQ potential: 24.79(20) MeV (GFMC) and 22.75(10) MeV (VMC). 

The first FY/AGS calculation that went beyond the truncated NN t-matrix approximation by fully 
considering the 3 s 1 - 3 c? 1 tensor force (NN potential of Yamaguchi form) was achieved by Fonseca [44] 
(see also Table 3). The results showed that a large part of the tensor force effect to the 4 He binding 
energy was already contained in a truncated t-matrix approximation. 

In the early phase of the 4N problem various ab initio methods were also applied to 4N scattering. 
We do not discuss the corresponding results in detail, but, at least, give a list of these early 4N scattering 
studies: [210] (RGM), [211] (HH), [212] (FY), [213] (FY), [214] (AGS-like), [215] (RGM). 

Here we just mention a few points. Tjon [213] made an interesting finding, namely that the effect of 
p-waves in the 3N subamplitudes is not negligible for low-energy 4N scattering. This was later confirmed 
in an AGS calculation [216], where it was shown that the effect is particularly strong in p- 3 He elastic 
cross sections, but also noticeable in dd — > p- 3 R and dd — > dd reactions. These studies of low-energy 
4N scattering reactions by Fonseca with the AGS approach were continued in Ref. [217], and, as we will 
see in Section 3.2.3, have been further pursued till today. 



3.2 4N Problem: The Modern Phase 

The more modern phase of 4N calculations started at the end of the 1980s, when not only the first 
GFMC calculations with realistic NN potential models emerged, but, in addition in the early 1990s 
quite a variety of different groups applying different methods appeared on the scene. Actually, the list 
of methods is quite long: (i) GEM [76, 218], (ii) DMC [219], (iii) CC [120], (iv) IDEA [220], (v) 
CHH [67], (vi) FY in configuration space [221], and (vii) FY in momentum space [33] (these methods 
are discussed or at least mentioned in Section. 2). All these various techniques were employed in 4 He 
bound-state calculations, whereas a new onset of 4N scattering-state calculations came somewhat later 
as will be discussed below. 
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Table 9: 4 He binding energy in MeV for the MTV potential with various methods. 



Method 


GEM [76, 218] 


IDEA [220] 


DMC [219]) 


FY [33] 


CHH [67] 




31.357 


30.98 


31.5 


31.36 


31.355 



3.2.1 The Four-Nucleon Bound State 

A first GFMC calculation with realistic NN potential models was carried out by Carlson [90] using the 
BonnR [89] potential and the Reid [187] potential reduced to V8 form. He obtained an underbinding of 
about 5 MeV, which is similar to the results of previous realistic calculations [119, 183]. In Ref. [90] it was 
also shown that plausible models for the 3NF can remove the underbinding. With time, improved models 
for NN potentials and 3NFs evolved. Since the GFMC technique cannot treat momentum dependent 
interactions, (see Section 2.5.1) the method remained restricted to r-space interaction models. Thus a 
GFMC 4 He binding energy result for a modern realistic force model exists only for the AVI 8 [96] NN 
potential. With an additional inclusion of the 3NF model UIX [139] a binding energy of 28.30(2) MeV 
was obtained [139], which is in agreement with the experimental result. 

In the momentum space FY calculation of Ref. [33] quite a variety of NN potential models were 
considered: local and non-local s-wave models with and without spin dependence, local central models, 
separable models with 3 s\- 3 d\ tensor force, 1 sq and 3 si- 3 g?i interactions from Reid potential. An extensive 
comparison with other results in the literature was made. In Table 9 we show the results with the MTV 
potential, where a benchmark had already been established before (see Table 8). Comparing with the 
benchmark one sees that most of the techniques lead to similar results. One could even say that the 
precision of the benchmark result was further increased (compare GEM, FY, and CHH results from 
Table 9 and ATMS result from Table 8). At this point it should also be mentioned that the convergence 
of the FY calculation is shown for an increasing number of channels and it is evident that excellent 
convergence is established (the maximal number of channels was 92, note that the MTV potential 
reduces the 4N problem to a four-boson problem with a considerably smaller number of channels than 
given in Table 1 for the four-fermion case). Of all the results listed in Table 9 the IDEA binding 
energy shows the strongest deviation, but IDEA is not a full ab initio method, since it can be seen 
as equivalent to a HH expansion limited to the lowest HH state (see Section 2.3.1). It should also be 
mentioned that the FY result for the MTI-III interaction agreed very well with the high-precision result 
from the coordinate space FY calculation of Ref. [221]. 

Kamada and Glockle extended their FY calculations to other potential models. In Ref. [34] the 1 Sq 
and 3 si- 3 di interactions from various realistic NN potentials were considered. In comparison to Ref. [33] 
the precision of the calculation was improved by increasing the number of channels considered (Table 2 
shows examples for some of the NN potentials). In addition it was shown that potentials that are fit 
to describe np or pp scattering data lead to somewhat different Tjon lines. A further improvement was 
made in Ref. [1] by taking into account NN forces in all NN partial waves with j < 3 (results given 
in the last column of Table 2). The resulting binding energies ranged from 23.45 to 28.11 MeV. It 
is interesting to note that the lowest value, which arises from the Reid potential, is close to the CC 
results from the 1970s mentioned in Section 3.1. The formalism to incorporate 3N forces into the 4N 
Yakubovsky scheme was laid out in Ref. [35]. Fully converged momentum space FY calculations with 
modern realistic forces and various 3NF models were performed in Refs. [39, 222]. The variation of the 
4 He binding energy with modern realistic NN forces was found to be considerably smaller than with 
the non-modern realistic models used in Ref. [1]. It was observed in Ref. [139] that by adjusting the 
3NF in order to describe the triton binding energy one is led to a 4 He binding energy rather close to 
the experimental value. This apparently diminishes the need for a 4NF. The smallness of the 4NF is 
also confirmed in Ref. [223], where the leading contribution of a chiral 4NF to the a-particle binding 
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Table 10: The 4 He expectation values (T) and (V) of kinetic and potential energies, the binding 
energy BE, and the root-mean-square (rms) radius with AV8' NN potential [139] calculated with various 
methods (from Ref. [226]). 



Method 


(T) [MeVl 

\ / L J 


(V) [MeV] 

\ / L J 


BE [MeV] 


J(r 2 ) [mil 


FY 


102.39(5) 


-128.33(10) 


25.94(5) 


1.485(3) 


GEM 


102.30 


-128.20 


25.90 


1.482 


SVM 


102.35 


-128.27 


25.92 


1.486 


HH 


102.44 


-128.34 


25.90(1) 


1.483 


GFMC 


102.3(1.0) 


-128.5(1.0) 


25.93(2) 


1.490(5) 


NCSM 


103.35 


-129.45 


25.80(20) 


1.485 


EIHH 


100.8(9) 


-126.7(9) 


25.944(10) 


1.486 



Table 11: Probabilities of 4 He total angular momentum components in % for AV8' potential calculated 
with various methods (from Ref. [226]). 



Method 


s-wave 


p-wave 


ci-wave 


FY 


85.71 


0.38 


13.91 


GEM 


85.73 


0.37 


13.90 


SVM 


85.72 


0.368 


13.91 


HH 


85.72 


0.369 


13.91 


NCSM 


86.73 


0.29 


12.98 


EIHH 


85.73(2) 


0.370(1) 


13.89(1) 



energy was estimated to be of the order of a few 100 keV. 

During the 1990s not only was great progress made with the GFMC and FY calculations, but also 
with some other of the above mentioned techniques (see e.g. CHH [71], GEM calculations for 4-body 
hypernuclei [224]). In addition other players entered the game with further methods: the SVM [78, 79], 
the NCSM [110, 225], and the EIHH methods [113, 114]. A manifestation of all this progress was 
the benchmark calculation of Ref. [226], where seven different theory groups applied seven different 
techniques to solve the 4 He bound-state problem with the realistic AV8' potential [139]. Besides the 
binding energy various calculational results were compared: expectation values for kinetic and potential 
energies, separate expectation values for various potential terms, the NN correlation function, and the 
probabilities of the total angular momentum components. In Tables 10 and 11 we list a few of these 
results. One sees that a very good agreement between the various methods was obtained. A closer 
inspection of Table 10 shows that deviations for the binding energy and rms radius are within 0.5%. 
Except for NCSM and EIHH similarly good agreement is found for (T) and (V). However, the NCSM 
results are still within 1% and EIHH within 1.5% of the others, but one should note that the EIHH results 
were calculated with the bare operators. The probabilities of the various wave function components 
displayed in Table 11 also show very good agreement apart from a small excursion in NCSM. Like the 
benchmark for the simple MTV potential from the early 1980s (see Table 8), the benchmark for the 
realistic AV8' of Ref. [226] can be considered as a milestone in the development of ab initio calculations 
for the 4N system. 

More recently, a similar agreement was achieved considering different NN potentials and 3NFs. As 
an example we show in Table 12 results for the interaction model AV18+UIX. 

A different approach for describing the nuclear force, namely without a 3NF, is made via a coupled 
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Table 12: 4 He binding energy in MeV for AV18+UIX Hamiltonian calculated with various methods. 



Method 


GFMC [227] 


FY [39] 


HH [228] 


EIHH [229] 




28.34(4) 


28.50 


28.46 


28.42 



channel calculation, where N and A degrees of freedom are treated on the same level by using an 
appropriate potential like the CD-Bonn+A potential [230]. Such a model, however, led to a slight 
underbinding in three- and four-nucleon systems implying that many-nucleon forces not accounted for 
by the A isobar make a non negligible contribution to the binding energies [231]. 

3.2.2 The Four-Nucleon Scattering State 

Continuum wave functions describing the 4N system in scattering states can be studied both in hadronic 
reactions and in reactions with external electroweak probes, i.e. photons, electrons and neutrinos. It is 
worthwhile pointing out that the information supplied by one type of reaction does not render redundant 
information gained from the other. Specific effects could have different relative sizes in the various types 
of reactions. In this context it is important that in electron scattering, as in principle also in neutrino 
scattering, energy transfer u and momentum transfer q can be varied independently. Thus, for a given 
u, i.e for a fixed continuum state, one can study the q dependence of various effects. In addition one 
should note that, due to the existence of two-body current operators, two completely phase equivalent 
potential models can nonetheless lead to different responses to external probes. 

3.2.3 Hadronic Reactions 

Starting in the late 1990s considerable progress was made in 4N continuum calculations. The zero-energy 
scattering of n- 3 H and p- 3 He was studied in a HH calculation by Viviani, Rosati, and Kievsky [232] 
through application of the Kohn variational principle (see Section 2.3). Together with the Coulomb force 
various versions of the Argonne potentials and various 3NF models were considered. The convergence 
of the results was checked by increasing the number of channels in the HH calculation. Almost simulta- 
neously Ciesielsky, Carbonell, and Gignoux began a study of the 4N scattering problem in a coordinate 
space FY calculation [233]. This enabled a successful comparison of HH and FY results for the n- 3 H 
singlet and triplet scattering lengths for the AV14 potential. Concerning other results in Ref. [232], it 
was found that the 3NF decrease the n- 3 H total zero-energy scattering cross section by about 10%, 
which then led to a good agreement with experimental data. A somewhat inferior agreement, with a 
difference of about 3%, was obtained for the coherent scattering length. 

The already above mentioned FY study [29] was an extension, by inclusion of spin and isospin 
degrees of freedom, of the four-boson configuration space calculation of Ref. [28] . The Grenoble group 
used the MTI-III NN potential without consideration of the Coulomb force in their first study. For 
energies below the three-body break-up threshold they calculated phase shifts for the elastic (3+1), 
as well as phase shifts for the inelastic (3+1)— >(2+2) reaction. A particularly interesting finding of 
their calculation was the good description of the resonant structure (p-wave resonance) of the nt total 
elastic cross section obtained with the simple MTI-III potential (see Fig. 8a). This fact initiated a very 
fruitful and interesting discussion further on. In an improved calculation [2] various 1 sq and 3 si- 3 d\ 
interactions from realistic NN potential models were used. In addition a phenomenological 3NF was fit 
to describe the binding energy of 3 H. Since the number of channels considered was not sufficiently high 
the calculation was reasonably well but not fully converged (a comparison of a corresponding calculation 
of the 4 He binding energy with precise ab initio results showed that less than 1 MeV was missing). We 
list the main findings of Ref. [2] : (i) only a rather weak dependence of the results on the choice of the 
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NN potential model, (ii) effect of the 3NF is very important at very low energies (as already found in 
Ref. [232], it led to a very good agreement with the zero-energy experimental cross section for elastic nt 
scattering), (iii) in contrast to the MTI-III case a good description of the cross section in the p-wave 
resonance region was not obtained with realistic NN potentials (see Fig. 8b), and (iv) 3NF effects are 
very small in the resonance region. 




Figure 8: n- 3 H total cross section: FY results [2, 29]. (a): for MTI-III, with 4N-s- (dotted), p- (dash- 
dotted), and (s +p)-wave (full); (b): with 4N-(s + p)-wave for AV14 (dashed) and AV14+3NF (full). 
Experimental data (crosses) [234]. 



The FY calculation for n- 3 H elastic scattering of Ref. [2] was further improved by Lazauskas and 
Carbonell [31] by including the NN force in NN partial waves with j < 3, and additionally also in 
3 / 4 - 3 /i 4 . As interactions both the AV18+UIX and the non-local NN potentials INOY [235] were used 
(the INOY non- localities were constructed such that an additional 3NF becomes obsolete). It was 
shown that the non-local model leads to rather good descriptions of the 4 He binding energy and of the 
zero-energy n- 3 H elastic scattering cross section, but fails, as does the AV18+UIX interaction model, 
to describe the cross section in the resonance region. 

At the end of the 1990s, the AGS approach for 4N scattering (see Section 3.1) was also improved [236]. 
The calculation was formally very similar to the previous work in Ref. [44] but contrary to the latter 
the interaction model became more realistic by inclusion of the interaction in NN p-waves. However, 
like the s-wave interactions the p-waves were also represented by separable expansions. In addition the 
Coulomb force was not included. Low energy elastic n- 3 H, (3 + 1) — > (2 + 2), and dd scattering were 
discussed. Here we only mention two results: (i) the calculated analyzing power A y in the 3 He(p,p) 3 He 
reaction was considerably smaller than the experimental A y) and (ii) the effect of the NN p-waves on 
the n- 3 H elastic scattering cross section was found to be different from Ref. [2]. 

In Ref. [87] Hofmann and Hale studied the 4 He compound system, i.e. the 4N system in the isospin 
T = channel, in the continuum at low energies. They performed both, a R-matrix analysis and a 
RGM calculation. In Ref. [237] their work is continued and applied to n- 3 H and p- 3 He scattering, where 
among other nuclear forces they also include the AV18+UIX model in the RGM calculation (interaction 
included for NN s-, p-, and partially also d- waves). Also in this case we want to mention two of the 
results: (i) the elastic n- 3 H cross section agrees with data in the very low-energy region, but as in 
Ref. [2] underestimates the resonance peak; and (ii) similar to Refs. [236, 238] the experimental A y 
in p- 3 He scattering is strongly underestimated (note that Refs. [237] and [238] were published almost 
simultaneously; further discussion of Ref. [238] is given below). 

As a matter of fact the n- 3 H elastic scattering results obtained with various methods (AGS, FY, 
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Figure 9: As Fig. 8 but for AV18 NN poten- 
tial calculated using three different meth- 
ods, almost identical results for FY and 
HH calculations (dotted) and AGS results 
(dashed); from Ref. [239]. 
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Figure 10: As Fig. 8 but calculated with various methods 
and potentials: AGS result [231] with CD-Bonn+A [230] 
(full), AGS result [240] with INOY04 [244] (dashed), FY 
result [239] with AV18+UIX (dotted). Data as in Fig. 8. 



HH, RGM) showed a large variance among each other for some observables (e.g. resonance cross section 
in n- 3 H scattering). Improvements in the AGS, FY, and HH calculations were made in order to produce 
a benchmark [239], wherein various scattering observables were calculated by the three participating 
groups. In Fig. 9 we show the cross section for elastic n- 3 H scattering as an example. It is evident 
that FY and HH results agree perfectly, whereas the AGS cross section is somewhat different. Similar 
findings could be made for a few others of the cases studied. It was then concluded that probably the 
AGS separable expansion of the NN potentials with just one separable term is not sufficient for p-waves. 

The above example of a benchmark test calculation is of course only possible for ab initio calculations, 
since only in this case does a unique ab initio result exist. The conclusion to be drawn from the 
benchmark case above was quite clear: improvements were necessary for the AGS calculation. In fact 
not much later the AGS calculation was considerably improved. In Ref. [240] Deltuva and Fonseca 
abandoned the various separable expansions usually made in the AGS formalism. Instead, they treated 
the NN t-matrix exactly and solved the resulting three-variable integral equations without any further 
reduction. In addition the number of NN, 3N, and 4N partial waves was increased up to convergence. 
The new calculation for low-energy elastic n- 3 H scattering confirmed in principle the results of Refs [2, 
29] concerning the resonance peak, although the cross section results were slightly higher, probably 
due to better convergence in the AGS calculation. Agreement was also obtained with the phase shift 
calculations of the Grenoble and Pisa groups given in the benchmark of Ref. [239]. 

A few years earlier the Pisa group extended the HH calculations for the zero-energy p - 3 He scatter- 
ing [232] to low-energy scattering employing the AV18+UIX nuclear force. Their results were published 
in Ref. [238] together with new accurate data of the proton analyzing power A y . It was found that 
the theoretical results for A y were up to 40% smaller than the data, thus leading to a situation which 
seemed to be similar to the well-known A y puzzle in 3N scattering. The 4N A y problem was confirmed 
in a later study [241], where new accurate experimental results were presented together with new HH 
calculations of the Pisa group. Shortly afterwards Deltuva and Fonseca came out with a fully converged 
AGS calculation for p- 3 He scattering [243]. Various NN potentials, including the non-local INOY po- 
tential, were employed. The inclusion of the Coulomb force was achieved by using a screened Coulomb 
force according to the methodology developed in [3]. Large Coulomb effects were found for all the 
studied observables. The calculation described experimental data quite well, except for the proton A y , 
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where the discrepancy between theory and experiment, found in [238, 241], was confirmed. 

The AGS calculations were extended to 4N scattering in the isospin T = channel [242], which 
makes the consideration of couplings between the n- 3 He, p- 3 H, and d-d channels necessary. The same 
NN potential models as in Ref. [243] were used. Observables pertaining to the six different elastic 
and transfer reactions were obtained and compared to experiment. The agreement between theory 
and experiment was fairly good, but for some observables stronger differences were present. In order 
to include many-body forces the study was continued using the CD-Bonn+A two-baryon model [231]. 
However, the inclusion of the A was unable to resolve the existing discrepancies with experimental data 
and the effect of the 4NF turned out to be very small. In Fig. 10 we show as an example the n- 3 H 
elastic scattering cross section with various interaction models. One sees that all of them fail to describe 
the p-wave resonance peak. 

In a further AGS calculation various observables of nucleon transfer reactions in low-energy d-d 
scattering were studied [245]. For example, it was found that dd fusion with spin-aligned deuterons is 
significantly suppressed at a few MeV but not in the keV regime (see Fig. 11). 

In continuation of their work Hofmann and Hale [88] carried out another, more complete, RGM 
calculation for the 4 He compound system using as nuclear force the model AV18+UIX. The resulting 
phase shifts were compared to a comprehensive i?-matrix analysis. The agreement between calculation 
and analysis was rather good in most cases and 3NF effects were found to be quite large for some specific 
cases. The calculation was restricted to a value of L = 3 for the maximal orbital angular momentum 
between the various two-body clusters, which, however, for some observables is not sufficient as pointed 
out later in Ref. [245]. 

Another RGM calculation was carried out by Quaglioni and Navratil [246] . They considered neutron 
scattering on 3 H (and also on 4 He and 10 Be, see Section 2.3.4) as well as proton scattering on 3 He (and 
also on 4 He). As opposed to the RGM calculations mentioned above their nuclear wave function 
consisted in a fragmentation of the form 3 X-N, with the 3N cluster in its ground state. The calculation 
is based on the NCSM, with corresponding HO expansions for the various wave function pieces to be 
determined. Several NN potentials were used, but a 3NF was not considered. The obtained phase 
shifts were compared to AGS results of the Lisbon group [240, 243] for the same interaction model. A 
fairly good agreement was found except for n- 3 H p- waves, where significant differences were observed 
at Ecm > 1-5 MeV. Certainly space for further improvement remained, which could be obtained by 
extending the RGM calculation beyond the single-fragment approximation. On the other hand we 
should mention that the 4N calculation merely served as a testing ground, because the principle aim 
expressed in Ref. [246] was the determination of phase shifts for elastic n- 4 He, n- 10 B, and p- 4 He 
scattering. 

In further calculations of the Pisa group various nuclear force models were considered, among them a 
chiral force (I-N3LO NN potential [128], N2LO-3NF [247, 248]). Their work consisted of high-precision 
calculations of the 4N bound state, the zero-energy scattering states [73] and of p- 3 He scattering [249]. 
In this study it was found that use of the chiral nuclear force gave improvements when comparing 
to experimental results. In fact a considerably better description of the analyzing power A y in p- 
3 He scattering was obtained. However, some differences between experimental and theoretical results 
remained as can been seen in Fig. 12, where we show A y results with conventional and chiral nuclear 
forces. One notes that there is a relatively strong increase of A y due to the chiral 3NF. 

The 4 He compound system exhibits a + resonance between the thresholds for p- 3 H and n- 3 He break- 
up. This resonance was studied by Lazauskas in a FY calculation for elastic p- 3 H scattering [32]. Various 
nuclear force models (Vnn + 3NF + Coulomb force) were used and effects from charge independence 
and charge symmetry breaking in the strong part of the nuclear interaction were taken into account. 
The convergence of the calculation was checked and a convergence level better than 0.5% was obtained. 
It was found that none of the tested Hamiltonians were able to reproduce simultaneously the shape of 
the resonance and the 4 He binding energy. 
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Recently another benchmark [40], a follow-up of that in Ref. [239], was performed for n- 3 H and 
p- 3 He scattering using AGS, HH, and FY techniques and taking various realistic NN potentials. The 
agreement among the three different approaches for the considered observables (phase shifts, differential 
cross sections, analyzing powers) was found to be better than 1%. Thus another milestone in the 
development of ab initio calculations in the 4N system was set. 

Elastic n- 3 H andp- 3 He scattering above the four-body break-up threshold was considered in a RGM- 
NCSM calculation limited, however, to a fragmentation into a nucleon and the 3N system in its ground 
state [252]. Nonetheless, a fairly good agreement with experimental results was obtained (differences 
were particularly present at forward angles). Only very recently a breakthrough was achieved by Deltuva 
and Fonseca [153] in an AGS calculation of the elastic n- 3 H scattering, considering the full four-body 
continuum. 

Finally we would like to mention that recently phase shifts of elastic N - 3 He scattering were calculated 
with realistic NN potential models for energies below the 3 He break-up threshold by applying bound- 
state like methods [253]. 




Figure 11: The quintet suppression factor 
(QSF) for d+d ->■ p+ 3 R and d+d ->■ n+ 3 He 
reactions as function of the deuteron energy 
in the laboratory system from an AGS cal- 
culation with the INOY04 potential [245]. 



Figure 12: A y of p- 3 He scattering at E p =4 MeV cal- 
culated with various interactions. HH results [250]: 
AV18+UIX (dashed), I-N3LO+N2LO (solid); AGS re- 
sults [243]: I-N3LO (dotted) and INOY (dash-dotted). 
Experimental data from Ref. [251]. 



3.2.4 Reactions with External Probes 

Among the various perturbation-induced reactions with 4 He, elastic electron scattering was naturally 
considered first in ab initio calculations. Since the total angular momentum of the 4 He ground state is 
zero, there exists only an elastic 4 He charge form factor. 

For inelastic perturbation-induced reactions one has to include the final state interaction (FSI) of the 
disintegrated 4N system. This can be accomplished either by a direct calculation of the 4N continuum 
state or by the methods described in Section 2.6. 

A calculation of the elastic 4 He charge form factor with a nuclear one-body charge operator was 
carried out by Tjon in a FY calculation in 1978 [188]. By taking for the NN interaction a truncated Reid 
potential in t 00 approximation (see Section 3.1) he obtained a rather good agreement with experimental 
data up to almost the first form factor minimum at about q 2 = 10 fm -2 . A more recent calculation 
was made in Ref. [254], where various modern nuclear Hamiltonians were considered and the 4 He wave 
function was computed via a HH expansion. Besides the one-body charge operator, two-body operator 



corrections (TBOC) were also taken into account. The dependence of the form factor on the force model 
was found to be rather weak. As shown in Fig. 13 there is a sizable TBOC contribution which then 
leads to good agreement with experimental data in the entire q range considered. This demonstrates 
that the modern nuclear force models describe the internal structure of 4 He with great precision. 




Figure 13: 4 He charge form factor: HH re- Figure 14: Cross section ratio 4 He(e, e'p) 3 H/ 
suits [254] with AV18+UIX potential; experimen- 4 He(e, e'n) 3 He at q=77 MeV/c. RGM results (see 
tal data from [255]. text) and experimental data from [267]. 

Of all the inelastic perturbative-induced reactions with 4 He, perhaps most attention was given to 
photodisintegration. Calculations were made for various low-energy break-up reactions, i.e. 4 He(7, p) 3 H, 
4 He(7,n) 3 He, and 4 He(7,cf) 2 H. For the latter case, there had already been early interest in the cross 
section at very low energies in order to determine the astrophysical S-factor. Because the reaction 
is of isoscalar nature, El transitions are strongly suppressed and the reaction is mainly driven by E2 
transitions. A review of the field up to 1988 is given by Weller and Lehman [256]. Later more ambitious 
calculations were carried out by applying the RGM technique (see Ref. [257] for a summary). With 
regard to the 4 He low-energy photodisintegration into the (3+1) channel, abundant experimental data 
are available although some of the cross section results are controversial. It is beyond the scope of the 
present review to describe the rather complicated situation, instead we refer to Ref. [258] (see also the 
discussion below). 

An important astrophysical reaction in the 4N system is the hep process (p + 3 He — > 4 He + e + + v e ). 
In this case the S'-factor was calculated by Marcucci, Schiavilla, and the Pisa group using HH expansions 
for bound and scattering states [259, 260]. A calculation of the thermal neutron capture on 3 He was 
pursued along similar lines [261]. In all these cases the nuclear force model AV18+UIX was used, but 
it should be pointed out that the electroweak operators were determined in two different ways, namely 
using a conventional nuclear physics approach [259, 261] and a chiral expansion [260, 261]. The following 
S-factors were obtained S hep (0) = 9.64 • 10~ 20 keVb [259] and Sh ep (0) = (8.6 ± 1.3) • 10~ 20 keVb [260]. 
The error estimate in the chiral calculation stems from the cut-off dependence inherent in the chiral 
approach (in Ref. [262] a doubling of the error due to the cut-off dependence is recommended). For 
the thermal neutron capture on 3 He the chiral expansion was carried out up to the 4th order (N 3 LO) 
and showed a rather surprising convergence pattern (strong effect of the N 3 LO term). A very good 
agreement with experimental data was attained. On the other hand one should mention that in order 
to determine the chiral Ml current operator four of the so-called LECs (low-energy constants) were 
fitted to electromagnetic nuclear observables (magnetic moments of 2,3 H, 3 He and np thermal capture 
cross section). As shown by Girlanda et al. [261] the conventional nuclear physics approach with a 
standard, though rather complicated, current operator, obviously without any open fit parameter, also 
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leads to a result rather close to the experimental data. 

Low-energy processes with 4 He were studied in electron scattering as well. In three experiments the 
+ resonance, which is located closely above threshold (see also Section 3.2.3) was investigated in the 
momentum transfer range 100 < q < 380 MeV/c [264, 265, 266]. By treatment of the resonance as a 
bound state a GEM calculation [263] with the AV8' NN potential and an ad hoc 3NF resulted in fairly 
good agreement with the experimental data. A further experimental study at low energy was made 
in Ref. [267] for the 4 He(e, e'p) 3 H and 4 He(e, e'n) 3 He reactions for excitation energies E x in the range 
22 < E x < 36 MeV and with q = 77 MeV/c. The experimental results were compared to those of an 
RGM calculation, which was carried out with a charge symmetric nuclear Hamiltonian with a semi- 
realistic Vnn consisting of central, spin-orbit, and tensor components. The Coulomb force was included 
and electromagnetic multipoles were considered up to the order 1 = 2. Even though the calculation did 
not consider three- and four-body break-up reactions (thresholds at E x = 26.1 and 28.3 MeV), a very 
good agreement of theoretical and experimental cross sections, both for proton and neutron knock-out, 
could be obtained. Because the uncertainty of the measured absolute cross section amounted to ±15%, 
it was safer to compare the cross section ratio R e = a(e, e'p)/cr(e, e'n). Of course very good agreement 
between theory and experiment was found for R e (see Fig. 14). This result for R e certainly did not 
leave much need for charge symmetry breaking terms in the nuclear Hamiltonian. This was in contrast 
to findings for the corresponding ratio R 1 obtained in 4 He photoabsorption, but, as already mentioned 
above, the 4 He photodisintegration cross sections were not yet sufficiently settled. 

Next we turn our attention to reactions beyond the low-energy regime. There, only very few theo- 
retical ab initio studies exist and in none of them does one find ab initio solutions of the 4N scattering 
problem beyond the three-body break-up threshold. It is true that continuum wave functions were 
calculated with AGS [268] and RGM [269] techniques up to rather high energies, but the many-body 
break-up was neglected in these calculations. Due to the recent success of Ref. [153] in treating the full 
four-body scattering, one can envisage in the not too far future a similar progress for the perturbation 
induced break-up of 4 He. For A > 4 the proper implementation of the full continuum is presently out 
of reach with ab initio methods, where the continuum wave function is calculated explicitly. Using, 
however, the approach discussed in Section 2.6.1 one reduces the continuum-state problem to a much 
simpler bound-state like problem, which then can be solved with standard bound-state techniques. 

A first use of integral transforms in few-nucleon ab initio calculations was made by Carlson and 
Schiavilla. They calculated the Laplace transform (Euclidean response) of the 4 He inclusive longitudinal 
response function Ri,{q, oo) via the GFMC technique [156] for various q- values. As nuclear forces the V8 
form of the AV14 potential and the UVIII-3NF [149] were employed. Comparisons with experimental 
results were made for the Euclidean response, which became possible by making Laplace transforms of 
the experimental Ri{q,oo) at a constant q and assuming in addition a specific high-energy behavior of 
Ri(q,u). For a more direct comparison of theory and experiment it is necessary to invert the integral 
transform. Since the GFMC calculation contains statistical errors, and as any numerical calculation 
also other small errors, the inversion of the Laplace transform is rather unstable. A similar instability 
for the inversion was found for the Stieltjes transform [151]. As described in Section 2.6.1, among the 
various proposed transforms the LIT turned out to be the most adequate one, since stable inversion 
results could be obtained. 

In the following years the LIT method was applied to quite a number of perturbation-induced reac- 
tions with three- and four-body nuclei [66] . Bound and LIT states were determined via HH expansions 
of CHH and/or EIHH type. The two very first LIT applications were made for 4 He. Almost simulta- 
neously calculations were carried out for the inclusive (e, e') response function Ri(q,u) [270] and for 
the 4 He total photoabsorption cross section [271] (see also Ref. [272]). The latter was calculated in un- 
retarded dipole approximation, which is very reliable in the low-energy region. In fact in [158] besides 
benchmarking FY and LIT calculations (shown in Fig. 7), a check of the validity of the unretarded 
dipole approximation was also made. It was found that retardation and higher multipole contributions 
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are quite small (below 2%) up to moderate photon energies (to < 50 MeV). 

The potential models used in both the virtual and real photon calculations were the rather simple 
MTI-III and similar central potential models. For Rl a fair agreement with experimental data was 
obtained in the considered momentum transfer range 300 < q < 500 MeV/c. It was found that 
FSI are quite important and that they are not negligible in the quasi-elastic peak, even at q = 500 
MeV/c [270, 273] (see Fig. 15, where in addition a result of a more realistic calculation with AV18+UIX 
is shown). 

The situation was quite different for the 4 He total photoabsorption cross section. As already men- 
tioned above the experimental results did not lead to a clear picture. In fact, the cross sections were 
rather different in the giant dipole resonance region. The LIT result agreed quite well with those data 
sets exhibiting a high cross section peak, while, on the other hand, data sets with a lower peak height 
were generally viewed more favorably in the 1980s and 1990s, with the exception, however, of results 
deduced from elastic photon scattering off 4 He [276]. Also due to the findings of the LIT ab initio calcu- 
lation new experimental interest arose and various new data became available with time [258, 277, 278]. 
But before coming to a comparison of theory with the new experimental results we would first like to 
add some further information. In the meantime LIT calculations for the 4 He total photon absorption 
cross section were also performed with modern realistic nuclear forces (AV18+UIX in Ref. [229], chi- 
ral 2NF and 3NF in Ref. [279]). EIHH [229] and NCSM [279] techniques were used. The results of 
both calculations were rather similar and confirmed the quite pronounced low-energy cross section peak 
observed in the first LIT calculation [271]. We should mention two further recent calculations of the 
4 He total photoabsorption cross section with similar results as in Refs. [229, 279], namely a LIT/EIHH 
calculation with a UCOM version of AV18 by Bacca [280] and a calculation by Horiuchi et al. [155] 
in a complex scaling approach with the SVM method employing the AV8' potential and an additional 
phenomenological 3NF. 

As pointed out in Section 2.6.1 the LIT method can also be applied to exclusive reactions. In 
fact total cross sections for the 4 He(7,n) 3 He and 4 He(7,p) 3 H [281] reactions were calculated with the 
MTI-III potential. A nice agreement was obtained for the sum of both cross sections (the only open 
channels at low energy) and the total low-energy photoabsorption cross section calculated with the 
same potential model and by using the LIT technique for inclusive reactions. At higher energies the 
difference of both cross sections, i.e. the cross section for three- and four-body break-up reactions, was 
in fair agreement with corresponding experimental data. 

Now we come back to the comparison of theory and experiment for the low-energy 4 He photodisin- 
tegration. In Fig. 16 we show the cross section of the reaction 4 He(7,p) 3 H from the LIT calculation by 
Quaglioni et al. [281] in comparison to the most recent experimental data. One sees very good agree- 
ment of the LIT results with the energy dependent cross section of Ref. [258], although the theoretical 
results overestimate the data by about 10%. However, there is a strong disagreement with the data of 
Ref. [278]. It should be mentioned that recently Ref. [258] has cautioned against the use of the majority 
of radiative-capture data, like those of Ref. [278] for 4 He photodisintegration. We want to remind the 
reader that the LIT calculation of Ref. [281] was made with the simple MTI-III potential. ^From a com- 
parison of the 4 He total photo absorption cross section calculated with the MTI-III potential [271, 272] 
and with modern realistic nuclear Hamiltonians [229, 279] one can estimate that the MTI-III result 
should be about 10% too high. Thus one may expect to find a rather good agreement of a modern 
realistic ab initio result with the cross section of Ref. [258] . 

Recently LIT calculations with the EIHH method were carried out [282, 283] for the inclusive 
response i?i(g,w) of 4 He with modern realistic nuclear Hamiltonians (AV18 + UIX-3NF, AV18 + TM'- 
3NF [284]). Note both Hamiltonians lead essentially to the same 4 He binding energy. A momentum 
transfer range between 50 and 500 MeV/c was considered. In particular, the 3NF effects on R L were 
investigated. For 300 MeV/c < q < 500 MeV/c 3NFs lead typically to effects not greater than 5- 
10%. For lower momentum transfers however much larger effects were found, e.g. at g=50 MeV/c 



r.P. 



-- TN 

■■■■ TN(PWIA) 
— AV18+UIX 
□ Bates 1988 
• Saclay 1994 




150 200 250 
co [MeV] 



350 



Figure 15: R L (q,oo) of 4 He at g=500 MeV/c: LIT/CHH 
result [270, 273] with semi-realistic TN potential (dashed), 
LIT/EIHH result [282] with AV18+UIX (full), PWIA re- 
sult (dotted) with TN potential in spectral function ap- 
proach (see Ref. [273]). Experimental data from [274] 
(open squares) and [275] (full dots). 
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Figure 16: 4 He(7,p) 3 H cross section: Ex- 
clusive LIT/CHH calculation [281] with 
MTI-III potential (full). Experimental 
data from [258] (full dots) and [278] (open 
diamonds). 



one observes reductions of more than 30% due to the 3NF. The size of the effect depends even quite 
significantly on the choice of the 3NF [283] (see Fig. 17). 

A first LIT calculation for the 4 He (e, e') transverse response function Rt{q, oj) was performed in 
Ref. [286]. The MTI-III potential was employed and a consistent meson exchange current (MEC) was 
constructed. In contrast to calculations of Carlson and Schiavilla who used the Laplace transform 
approach with more realistic nuclear forces [287, 288], no strong MEC effect could be found with the 
semi-realistic MEC of the MTI-III potential at q = 300 MeV/c. 

In addition to the inclusive (e, e') reactions off 4 He the following two exclusive reactions were cal- 
culated using the MTI-III potential and CHH expansions: 4 He(e, e'j9) 3 H [160] and 4 He(e, e'ci) 2 H [161]. 
Last but not least we would like to mention a LIT benchmark calculation performed with NCSM and 
EIHH techniques [289], where the test was made for various responses employing a simple NN potential. 

Before concluding this section we have to mention that ab initio calculations were also performed 
for inelastic neutrino scattering off 4 He. These reactions play a crucial role in supernova explosion 
scenarios. In the LIT calculation of Gazit and Barnea [290] the interaction models AV8', AV18, and 
AV18+UIX were used and HH expansions were employed with the EIHH technique. For the electroweak 
part one-body operators were taken into account, but such that vector MEC were included via the so- 
called Siegert operators. As a final result the temperature averaged cross section was given as function 
of the neutrino temperature. 



3.3 Energy Spectra of A > 4 Nuclei 

Another milestone in the development of ab initio calculations in few-body physics was set in a series of 
articles by the Urbana-Argonne-(Los Alamos) collaboration who by using the GFMC approach paved 
the way for bound-state calculations of A > 4 nuclear systems. The first study was performed by 
Pudliner, Pandharipande, Carlson, and Wiringa in 1995 [291]. As interactions they basically chose 
the AV18 NN potential and the UIX-3NF. However the V8 part of AV18 together with the 3NF were 
calculated with GFMC while the remaining V9 - V14 and Coulomb parts were considered in first order 
perturbation theory. Further electromagnetic corrections were taken from VMC. 




Figure 17: R L {q,oo) of 4 He at various values of q: LIT/EIHH results [282, 283] with AV18 (dot-dashed), 
AV18+TM' (dashed), and AV18+UIX (full) potentials; experimental data from Ref. [285]. 



Table 13: Experimental and GFMC (AV18+UIX) energies of A = 6, 7 nuclei in MeV; also given is A, 
the deviation in MeV of the GFMC energy from the experimental one (from Ref. [139]). 



A X(J n ;T) 


GFMC 


Expt. 


A 


A X(J*-T) 


GFMC 


Expt. 


A 


6 He(0+; 1) 


-27.64(14) 


-29.27 


1.63 


7 He(|-;|) 


-25.16(16) 


-28.82 


3.66 


6 He(2+; 1) 


-25.84(11) 


-27.47 


1.63 


7 Li(F;|) 


-37.44(28) 


-39.24 


1.80 


6 Li(l+;0) 


-31.25(11) 


-31.99 


0.74 




-36.68(30) 


-38.76 


2.08 


6 Li(3+;0) 


-28.53(32) 


-29.80 


1.27 




-31.72(30) 


-34.61 


2.89 


6 Li(0+; 1) 


-27.31(15) 


-28.43 


1.12 


Ml"; 5) 


-30.88(35) 


-32.56 


1.68 


6 Li(2+;0) 


-26.82(35) 


-27.68 


0.86 


7 Li(r;|) 


-24.79(18) 


-28.00 


3.21 


6 Be(0 + ;l) 


-25.52(11) 


-26.92 


1.40 











Rather good results for the binding energies of the A = 6 nuclei were obtained, namely 28.2(8) 
MeV ( 6 He) and 32.4(9) MeV ( 6 Li), to be compared with the experimental values of 29.3 and 32.0 MeV, 
respectively. Also the calculated 6 Li proton point radius of 2.41 fm agreed well with the experimental 
value of 2.43 fm. In addition, the 1/2 - and 3/2~ resonances of 5 He were investigated, however, the 
calculated splitting of 0.8(3) MeV was much smaller than the observed 1.4 MeV. 

The GFMC study was continued in Ref. [139]. The calculation was extended to A = 7. The 
same interaction model (AV18+UIX) was used and quite a number of results were presented. Here we 
illustrate only the excitation spectra of the A = 6, 7 nuclei (see Table 13). It is evident that the energies 
are close, but somewhat above the experimental numbers. The authors believed that the discrepancy 
was most likely due to the phenomenological short-range part of the 3NF. In the next GFMC study 
also the A = 8 nuclei were included [227], although the nuclear force was still described by the model 
AV18+UIX. However now due to a constrained-path algorithm for the complex spatial and spin-isospin 
wave functions the precision of the calculation was improved. The authors estimated that, for a given 
Hamiltonian, binding energies could be obtained with an accuracy of 1 to 2 %. In Table 14 we list the 
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Table 14: Experimental and GFMC (AV18+UIX) energies of nuclear ground states in MeV; also given 
is A, defined as in Table 13 (from Ref. [227]). 



A X(J^;T) 


6 He(0+; 1) 


6 Li(l+;0) 


7 He(§-;i) 


7 Li(l~;|) 


8 He(0+;2) 


8 Li(2+; 1) 


8 Be(0+;0) 


GFMC 


-28.11(9) 


-31.15(11) 


-25.79(16) 


-37.78(14) 


-27.16(16) 


-38.01 


-54.44 


Expt. 


-29.27 


-31.99 


-28.82 


-39.24 


-31.41 


-41.28 


-56.50 


A 


1.16 


0.84 


3.03 


1.46 


4.25 


3.27 


2.06 



Table 15: Experimental energies of A = 9, 10 nuclei in MeV; also given is A(3NF-model), the deviation 
in MeV of the GFMC energy calculated with AV18 and a specific 3NF-model from the experimental 
energy (from Ref. [292]). 



A X(J W ) 


Expt. 


A(No-3NF) 


A(UIX) 


A(IL2) 


A(IL3) 


A(IL4) 


9Li (D 


-45.34 


11.64(30) 


4.44(30) 


-0.66(40) 


-1.36(50) 


-2.26(40) 


9 Li(D 


-42.65 


8.65(30) 


3.25(30) 


-1.15(40) 


-1.55(40) 


-1.85(50) 


9Li (D 


-39.96 


7.86(30) 


2.06(40) 


-1.14(40) 


-1.04(40) 


-0.64(40) 


9 Li(F) 


-38.91 


9.21(30) 


3.71(30) 


-0.09(40) 


0.21(40) 


-1.89(40) 


9 Be(f) 


-58.16 


11.76(40) 


3.06(30) 


-0.04(50) 


0.36(50) 


0.16(60) 


9 Be(|-) 


-55.73 


12.23(30) 


4.43(50) 


-0.07(50) 


1.63(40) 


0.63(50) 


9 Be(f-) 
10 Be(0+) 


-55.36 


10.36(40) 


4.46(60) 


1.06(40) 


-0.34(50) 


-0.14(40) 


-64.98 


12.98(50) 


5.78(60) 


-1.82(70) 


-0.22(70) 


-2.42(60) 


10 Be(2+) 


-61.61 


13.91(50) 


4.51(60) 


-0.19(50) 


1.71(60) 


0.51(60) 


10 B(3+) 


-64.75 


16.15(60) 


5.75(40) 


-0.85(50) 


0.65(50) 


-0.85(60) 


10 B(1+) 


-64.03 


12.43(60) 


3.73(50) 


-0.67(40) 


1.23(50) 


-1.07(50) 



calculated ground-state energies of the various nuclei considered. Comparing the results for the A = 6, 7 
nuclei with those of Table 13 one finds for three of the four cases that the theoretical results moved 
closer to the experimental values, however discrepancies still remained. In Ref. [227] a large number of 
additional results are given, among them being the spectra of various nuclei. A particularly interesting 
finding was the strong 2a clustering in the 8 Be ground state. In the conclusion the authors stated that 
the general features of the light p-shell nuclei were described fairly well, including the binding energies, 
and the correct ordering and approximate spacing of the spectra. They further noted specific failures 
of the AV18+UIX Hamiltonian, namely a few percent underbinding in N = Z nuclei (which increases 
as the n-p asymmetry increases) and insufficient spin-orbit splittings. 

In order to improve the comparison with experiment new 3NFs were developed, namely the five 
Illinois models (IL1-IL5) [293]. Compared to the UIX-3NF, which contains two terms, the Illinois-3NFs 
include two additional terms, a two-pion exchange s-wave term and three-pion ring term. The overall 
strengths of the four terms together with a cut-off factor in the radial dependence, were adjusted to 
obtain best fits to the energies of 17 narrow states in 3 < A < 8 nuclei (at most three parameters at a 
time could be uniquely determined and thus various models were constructed). The effect of the Illinois 
3NFs was then studied in a GFMC calculation for the A = 9, 10 nuclei [292]. In Table 15 we show 
some of the results obtained. Without a 3NF one observes a relatively large underbinding of about 8 to 
16 MeV which reduces to an underbinding of about 2 to 6 MeV when the UIX-3NF is introduced. The 
Illinois 3NFs, however, lead to much better results. The overall best result, including the 17 narrow 
states, was obtained with the IL2-3NF. In fact the following root mean square deviations from the 
experimental energies were found: 0.6 MeV (AV18+IL2), 0.8 MeV (AV18+IL3), and 1.0 (AV18+IL4). 



Table 16: 6 Li binding energy for Volkov and Minnesota [294] (MN) potentials with various ab initio 
methods (note MN results are taken from Ref. [298]) 



Potential 


SVM [78] 


HH [61] 


HH [64] 


NCSM [298] 


EIHH [113] 


Volkov 


66.25 


66.57 


66.49 






MN 


36.51 






36.50(4) 


36.64(7) 



This contrasts with rms deviations of 3.5 MeV for AV18+UIX and 9.9 MeV for AV18 alone. A further 
GFMC calculation using the AV18+IL2 Hamiltonian studied states with higher excitations in A = 6 - 8 
nuclei [140]. This calculation benefited from a significant improvement of the VMC trial functions. A 
good representation of the experimental spectra was found i.e. the rms deviations from 36 experimental 
energies with 6 < A < 8 was only 0.6 MeV. In a later calculation even the 12 C ground-state energy was 
calculated with the AV18+IL2 Hamiltonian [4]. Again, use of only the AV18 resulted in an underbinding 
of about 20 MeV which was reduced to about 3 MeV with the addition of the IL2-3NF. 

Other ab initio methods were not on a comparable level with the GFMC calculations for quite some 
time. The SVM technique was applied to A = 2 - 7 nuclei [78] using simple NN potentials, like the MTV 
and Volkov models. Similar simple potential models were taken in a HH calculation [61] based on an 
improvement of the IDEA approach (see Section 2.3.1). In Table 16 we list the 6 Li binding energy 
for two of these simple potential models as well as other available results from later ab initio studies. 
Comparing the various results of Table 16 one finds a rather good agreement. 
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Figure 18: Convergence pattern of 6 Li-NCSM 
ground state energy (squares) with AV8' and 
Minnesota [294] (MN) potentials (from [298]), 
also given results for GFMC/AV8' (full), 
SVM/MN (dotted), EIHH/MN (dashed). 



Figure 19: As in Fig. 18, but for 6 Li, 8 Be, 
10 B. NCSM results (AV8') with 2- (full circles) 
and 3-body (open squares) effective interaction 
(from [107]); also given GFMC results (full lines) 
with statistical error. 



It is interesting to investigate further progress in A > 4 calculations with the various ab initio 
methods of Table 16. The SVM technique and the HH approach of Ref. [61] were not developed further 
to treat the more realistic potential models. The HH method of Ref. [64] is quite recent and one has to 
await further developments. The EIHH approach was applied in several studies with central potential 
models in A = 6, 7 nuclei (see Section 3.4). In addition more realistic potentials were considered, e.g. 
6 Li with AV8' [115], although the results were not fully converged, particularly in the presence of a 
3NF. The approach which has seen the greatest progress is undoubtedly the NCSM, first presented by 
Navratil, Vary and Barrett [5] and recently reviewed in [112]. 



Table 17: Comparison of GFMC and NCSM energies of various nuclei calculated with the AV8' 
potential; listed in MeV are the GFMC energies and A = £(NCSM)-£(GFMC). NCSM results are 
calculated with two-body and partially also with three-body effective interactions, V2- e & and V3_ e g, 
respectively; number of considered HO quanta iV max in square brackets (from Refs. [107, 296, 299]). 



A X(J^;T) 


GFMC 


A(y 3 -eff) 


A(V 2 _ cff ) 


A X(J") 


GFMC 


A(y 3 -eff) 


A(y 2 _ efr ) 


e Li(l+;0) 
6 Li(3+;0) 
6 Li(0+;l) 


-28.19(5) 
-24.98(5) 
-24.15(4) 


-0.42 [6] 
-0.56 [6] 
-0.61 [6] 


-0.41 [10] 
-0.60 [10] 
-1.05 [10] 


8 Be(l+;0) 
8 Be(3+;0) 
8 Be(2+; 1) 


-32.77(15) 
-31.23(15) 
-32.7(1) 


1.64 [4] 

1.65 [4] 
0.86 [4] 


-0.03 [8] 
0.24 [8] 
-0.70 [8] 


6 Li(2+;0) 
8 Be(0+;0) 
8 Be(2+;0) 
8 Be(4+;0) 


-24.12(4) 
-47.89(11) 
-45.62(11) 
-38.69(11) 


0.15 [6] 
-0.57 [4] 
0.82 [4] 
2.64 [4] 


-0.11 [10] 
-1.83 [8] 
-0.47 [8] 
1.44 [8] 


9 Be(|-;I) 
10 B(l+;0) 
10 B(3+;0) 


-49.9(2) 
-55.67(26) 
-53.23(26) 


-0.52 [4] 
-1.60 [4] 


-0.30 [8] 
0.30 [8] 
-0.60 [8] 



In various NCSM studies [107, 295, 296, 297, 298, 299] spectra of A > 4 nuclei were calculated 
for various potential models, among them also the AV8' potential enabling comparisons with GFMC 
results of the Urbana-Argonne-(Los Alamos) group (see Table 17). Inspection of Table 17 shows that 
the NCSM results are rather close to those of the GFMC calculation. Thus, one may consider the 
NCSM results as almost, but not yet fully, converged. To investigate this better, we show in Fig. 18 the 
NCSM convergence pattern of the 6 Li ground-state energy with a two-body effective interaction. One 
sees that, in order to have a convergent result, the number of considered HO quanta, iV max , has to be 
quite large and certainly larger than the values given in Table 17. As shown in Fig. 19 the convergence 
can be accelerated if a three-body effective interaction is used instead of a two-body one. In fact, 
by inspecting the results of Table 17 one finds that to achieve the same quality as given by the pure 
two-body effective interaction result in general a lower iV max is needed when the three-body effective 
interaction is included. 

As evident from Fig. 18 the NCSM calculations were also performed with NN potentials other 
than the AV8'. The initially preferred model was the CD-Bonn potential, since 3NF effects are less 
pronounced in this case i.e. neglect of the 3NF has less important consequences. Later on, the chiral 
I-N3LO potential [128] was used in addition [298]. The first consideration of a 3NF in the NCSM 
calculations for A > 4 nuclei was made in Ref. [297], where the interaction model AV8'+TM' was used 
and accordingly a three-body effective interaction was employed. Taking A max = 6 for A = 6, 7 nuclei 
and iV max = 4 for A > 7 nuclei, NCSM calculations were carried out up to 13 C. Though such a relatively 
low value of iV max did not lead to convergent results (see also comparison of GFMC and NCSM results 
with AV8'+TM' in Ref. [300]), interesting 3NF effects like e.g. an enhancement of spin-orbit effects were 
found. In the large basis NSCM calculation (iV max = 8, 9) of 9,11 Be [299] the non-local NN potential 
INOY, which makes a 3NF obsolete, was added to the list of potentials. The INOY results suggested 
that a realistic 3NF could have an important influence on the parity inversion of the n Be ground state 
( n Be has an unnatural l/2 + ground state). 

In conjunction with the I-N3LO NN potential, two chiral 3NFs, models 3NF-A and 3NF-B, were 
used in the NCSM calculation of Ref. [111]. The excitation spectrum of 7 Li was calculated by taking 
Amax — 6, which did not allow fully convergent results to be reached. For the 7 Li binding energy the 
following results were obtained, 38.0 MeV (3NF-A) and 36.7 MeV (3NF-B), whereas one has a GFMC 
result (model AV18+IL2) of 38.9 MeV and an experimental value of 39.2 MeV. Thus the chiral forces, 
though having a soft core, did not lead to overbinding. As the authors pointed out, it was previously 
believed that only the addition of a repulsive core in the 3NF, as in the Urbana and Illinois 3NF models, 
could prevent overbinding. In addition, it was confirmed that the spectra depend on the choice of the 



specific 3NF. The study of NCSM calculations with chiral forces was continued in Ref. [301]. The chiral 
3NF was somewhat different from Ref. [Ill], in that the 3NF was chosen to be local and the constraint 
to fit the 4 He binding energy was loosened somewhat by considering additional observables of A > 4 
nuclei in order to determine one of the 3NF parameters. Besides these other observables, the spectra of 
A = 10- 13 nuclei were calculated using iV max = 6. Overall, the 3NF contributed significantly to improve 
the agreement between theory and experiment. This was especially well demonstrated in the odd mass 
nuclei for the lowest few excited states. 

Next we report on a recent NCSM study [302] which employed an importance truncation scheme, 
where a criterion based on perturbation theory permitted certain states to be discarded. In this way a 
NCSM calculation with iV max = 12 became feasible even for 12 C and 16 0. In addition a SRG-transformed 
chiral Hamiltonian was employed (I-N3LO and a N2LO-3NF version, where the 3NF parameters are 
constrained by both the triton binding energy and the triton /3-decay [303]). As described in Section 2.4 
the SRG introduces a prediagonalization of the Hamiltonian. Since it is important to understand the 
SRG transformation correctly, we reiterate here the essential features. In order to have a unitary 
transformation at the A-body level one should consider A-bodj operators in the SRG transformation. 
However, in practice the SRG transformation is performed only at the two- or at maximum the three- 
body level. A SRG transformation at the two-body level leads to a phase shift equivalent NN interaction. 
Nonetheless, depending on the bare interaction and on the particle system under consideration, one 
may obtain a unitary transformation at the two-body (three-body) level, which can lead to results 
different from those of the untransformed Hamiltonian for A > 2 (A > 3) systems. The unitarity of 
the transformation on a specific N-body level can be checked by the independence of the results on 
the flow parameter a. In other words, if one applies a SRG transformation at the two-body (three- 
body) level and obtains an a dependent binding energy in a calculation of an A-body system, then the 
transformation is not unitary, but induces a dependent three-body (four-body) and/or higher order 
many-body forces. Note that the case of a independence does not necessarily mean one is working with 
the original Hamiltonian, but that the induced many-body forces are a independent. On the other hand 
one is certain to be working with a NN phase shift equivalent interaction. It should be obvious that in 
presence of a bare 3NF the SRG transformation should be made at least at the three-body level. 

We come back now to the calculation of Ref. [302] where in fact the a dependence of the results was 
checked. By employing only the NN potential and the SRG transformation at the two-body level, it 
was shown that the ground-state energies of the nuclei considered ( 4 He, 6 Li, 12 C, 16 0) had a strong a 
dependence thereby indicating that the transformation was not unitary on the A-body level. By bringing 
the SRG transformation to the three-body level the results became approximately a independent at 
iV max = 12. It is interesting to note that an 16 binding energy of 119.7(64) MeV was obtained which is 
close to recent results of 121.0 and 119.5 MeV with CC [95] and UMOA [108] techniques, respectively 
(for the latter see the end of Section 2.4.1). With the additional inclusion of the bare 3NF it was found 
that only the 4 He binding energy is a independent. However by extrapolating to larger iV max it was 
also possible to determine the binding energy of 6 Li. On the contrary, even with extrapolations, there 
still remained a pronounced a dependence for 12 C and 16 O showing that the induced a dependent 4N 
and higher-order forces are quite large. 

The use of soft interactions makes full configuration NCSM calculations feasible. An example is found 
in Ref. [304], where A = 6, 8, 12 and 16 binding energies (and a few resonant states) were obtained. This 
was possible using the JISP potential [13] and an extrapolation procedure that was shown to be reliable 
in A = 2, 3, 4 test cases, where fully converged results can be obtained. Here one should also mention an 
interesting and promising exploration presented in Ref. [305] regarding the applicability of the MCSMD 
(see Section 2.5.3) to A > 4 systems. There a benchmark of the binding energies, point-particle rms 
matter radii, and electromagnetic moments for light nuclei ranging from 4 He to 12 C showed that results 
from full configuration interaction and MCSMD agree within a few percent at most. 

In recent years great progress was also made in the CC approach. In fact just mentioned above 
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was the CC result for 16 with the I-N3L0 potential. We already reported in Section 2.4.3 that CC 
binding energy results of 4 He were successfully benchmarked with FY results. The CC method is better 
suited for treating closed shell and neighboring nuclei and can be applied to nuclei with considerably 
larger A than with other ab initio methods. In Ref. [95], besides the above mentioned 16 case, even 
40 Ca and 48 Ca were studied with the I-N3LO NN potential in a A-CCSD(T) calculation. The model 
space comprising of 18 major HO shells was sufficiently large. For 48 Ca a similar underbinding per 
nucleon (0.4 MeV) as for 16 was found, whereas for 40 Ca a small overbinding of 0.08 MeV per nucleon 
was obtained. Thus a rather different effect of the 3NF would be needed for these nuclei. Further CC 
applications were performed mainly for A > 16 nuclei (see e.g. [306]), which is certainly out of the 
scope of the present review. First attempts of CC calculations in open shell nuclei with two nucleons 
(two holes) in the last shell were pursued by considering the 6 He nucleus [307] and working in small 
model spaces. The results were quite encouraging and prospects are good for future calculations in 
larger model spaces. 

Also results from lattice simulations became available lately (see Section 2.5.2). A LCEFT calcula- 
tion of the spectra of the three lightest a-nuclei was carried out in Ref. [308]. Particularly interesting 
is the result for the Hoyle state, where an energy of -85(3) MeV was found, which is close to the 
experimental value of -84.51 MeV. 

Another very promising approach for the ab initio treatment of a larger number of particles is the 
AFDMC technique mentioned in Section 2.5.1. For finite nuclear systems, however, applications have 
only be made with simpler interaction models like the AV6' NN potential [309]. 

3.4 Further Observables for A > 4 Nuclei 

In addition to excitation spectra further observables exist for nuclear bound states. One has charge 
and matter radii, magnetic and quadrupole moments, and electroweak transition strengths. Here we 
do not discuss these observables in detail but just comment on a few selected points. Since in most of 
the ab initio calculations wave functions of the various nuclear states are calculated, it is in principle 
straight forward to evaluate the above mentioned observables. In fact in most of the NCSM calculations 
of the previous section such observables were calculated. In addition we mention here two further 
NCSM studies: Ref. [310] (charge radii of 4 ' 6,8 He) and Ref. [311] (charge radii and electromagnetic 
moments of Li and Be isotopes). In GFMC calculations no wave functions are obtained, but nuclear 
radii and electromagnetic moments can be calculated nonetheless (see e.g. Ref. [227]). The GFMC 
calculation of electroweak transition matrix elements is described in Ref. [312], where quite a number 
of transitions were calculated with the AV18+IL2 interaction model. In most cases a good agreement 
with experiment could be obtained. The calculation was further improved by including MEC [313] when 
magnetic moments of A = 2 - 7 ground states and Ml transitions in A = 6, 7 nuclei were calculated. The 
contributions from MEC were considerable and led to an overall good agreement with experiment. 
Relevant MEC effects were also found in a HH hybrid calculation [314] for the 6 He /3-decay using 
the JISP potential and constructing a chiral weak current. The obtained decay rate was close to the 
experimental value (3% difference). 

Recently the radii of He isotopes were determined experimentally [315, 316] and compared to results 
from ab initio methods. In a correlation plot of the 6 He charge radius versus the two-neutron separation 
energy only the GFMC calculation (the only one with an explicit 3NF) agreed well with experiment. 
In a further rather recent experiment excitation energies of two 2 + excited states of 10 Be as well as 
the corresponding B(E2) values to the ground state were determined [317]. These data were compared 
with the results of a GFMC calculation that employed the AV18+IL7 interaction model (the IL7-3NF 
has to be understood as an improved IL2-3NF). A rather good agreement was found for three of the 
four observables, while the B(E2) value from the higher of the two excited states turned out to be 
considerably smaller than experiment. 



Now we turn our attention to reactions of nuclear systems with A > 4. First we consider n-a 
scattering which was calculated in Ref. [142] in a GFMC calculation with the interaction models AV18, 
AV18+UIX, and AV18+IL2. We depict those results in Fig. 20, where the calculated phase shifts of 
the 2 S\/2 and 2 p%/2 n-a partial waves together with a i?-matrix fit are shown. For the s-wave phase 
shifts one already has a good description with AV18 alone, whereas this is not the case for the p-wave. 
For the latter a very interesting 3NF effect was found. An inclusion of the UIX-3NF did not change the 
AV18 results much, whereas, at least for low energies, the phase shift moved very close to the i?-matrix 
fit when the IL2-3NF was taken into account. The nice success of a GFMC calculation for a low-energy 
scattering reaction in principle opened the door for further low-energy scattering calculations of other 
nuclear systems and of corresponding calculations of electroweak cross sections. 
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Figure 20: n- 4 He 2 Si/ 2 and 2 p3/2 phase shifts: GFMC 
results with AVI 8+UIX (triangles) and AV18+IL2 
(squares) from [142]; RGM/NCSM results with I- 
N3LO (full curves) from [318]. Plus signs: R- matrix 
analysis results from [319]. 



Figure 21: n- 4 He differential cross sec- 
tion (left panel) and analyzing power (right 
panel) at £^ b =17 MeV. RGM/NCSM results 
(full curves) from [320]; experimental data 
from [321]. 



In Fig. 20 we also show the I-N3LO results from a RGM calculation with the NCSM approach by 
Quaglioni and Navratil [318] who also used other NN potentials (CD-Bonn, V[ ow fc)- The potential model 
dependence was found to be rather mild. Only the 2p 3 / 2 phase shift computed with V\ ow k deviated a 
bit from the results arising from the two other potentials. Fig. 20 shows good agreement with the 
i?-matrix fit for the s-wave, whereas the 2p 3 / 2 phase shifts are too small, which as a comparison with 
the GFMC results show, is because a proper 3NF is missing. The RGM calculations were also carried 
out for p-a scattering. Similar comparisons with the i?-matrix fits were obtained i.e. agreement for 
the s-wave and disagreement for the p-wave. We should mention that the NCSM expansion included 
states up to iV max = 17, which provided nice convergence. Besides the 4 He ground-state, excitations of 
the a particle were also considered. No visible effect was found for the s-wave, while larger effects were 
present for the p-wave. In Ref. [318] further calculations were performed for the s-wave phase shift of 
n- 10 B scattering. Due to lack of data no comparison to experiment was possible but on the other hand 
the model space was not sufficiently large (iV max = 7) to reach convergence. 

The convergence was improved in a following NCSM study by Navratil, Roth, and Quaglioni [320]. 
They used the importance truncated NCSM mentioned in Section 3.3 and a SRG evolved (on the two- 
body level) Hamiltonian (I-N3LO). In this way the interaction was soft enough to reach convergence 
with iV max = 14-16 and thus nucleon scattering on various nuclei ( 4 He, 7 Li, 7 Be, 12 C, 16 0) could be 
investigated. The n-a scattering was used as a test case. In brief they found that (i) very good 
convergence was obtained for the phase shifts (in the iV max range from 11 to 17 the results looked 
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almost identical) and (ii) the calculation was benchmarked with the full NCSM and a nice agreement 
was obtained. The comparison to the i?-matrix fit did not change much, with agreement obtained for 
the 2 s 1 / 2 , 2 Pi/2, 2 dz/2 cases and, because of a missing 3NF, disagreement for the 2 p^/2 case. Comparison 
to experiment was made at E n — 17- 19 MeV, where 3NFs become less important. In Fig. 21 we show 
the results at E n = 17 MeV where it is seen that the differential cross section agrees very well with data. 
An overall good agreement was obtained for A y , but some differences were also present. Comparisons 
to data with very similar outcomes could be made for p-a scattering as well. In the case of n- 7 Li 
scattering a successful benchmark with the full NCSM was made again. Besides the 7 Li ( 7 Be) ground 
state two excited states (l/2~, 7/2 + ) were also considered in the RGM calculation. Various phase shifts 
for the (7+1) scattering were calculated, and for the reaction 7 Li(n, n') 7 Li(l/2~) even a comparison 
with total cross section data could be made at very low energy and a fair agreement was obtained (see 
Fig. 22). With an additional 3NF the comparison to data could probably be improved further. In the 
case of N- 12 C and N- 16 scattering, benchmarks with the full NCSM were no longer possible, since 
only the importance truncated NCSM could reach sufficiently high A max values. The SRG parameter 
A = a -1 / 4 was increased from 2.02 fm _1 to 2.66 fm _1 , the higher value leading to a shorter evolution and 
less soft potentials. The use of a small A results in large overbinding of heavier nuclei and a significant 
underestimation of their radii. We remind the reader that the transformed potential remains phase 
shift equivalent with the bare interaction, but that a dependent three- and higher many-body forces 
are induced (see discussion of Ref. [302] in Section 3.3). It is interesting to list the differences of the 
binding energies obtained with such a SRG evolved Hamiltonian and the corresponding experimental 
values. Underbinding of 0.3, 1.04, and 8.7(8) MeV was found for 3 H, 4 He, and 12 C respectively while 
overbinding of 11.4(8) MeV was found for 16 O. Here we discuss neither further calculational details nor 
the results for the obtained phase shifts of n- 12 C, p- 12 C, and p- 16 0, but refer the interested reader 
directly to Ref. [320]. We only mention that more difficulties were encountered in the case of 16 and 
that the authors conclude that the inclusion of additional excited states of the target nucleus would be 
beneficial in all studied systems particularly as A increases. 




Figure 22: Inelastic 7 Li(n,n') 7 Li(l/2-) 
cross section from RGM/NCSM calcula- 
tion of Ref. [320] with SRG/I-N3LO NN 
potential; SRG parameter A=2.02 fm _1 . 
Experimental data from Ref. [322]. 



Figure 23: 3 H(rf,n) 4 He S-factor from NCSM/RGM calcu- 
lation of Ref. [325] with SRG/I-N3LO NN potential in- 
cluding various numbers of deuteron virtual excitations 
in 3 si- 3 di (d*) and 3 d 2 (d'*) channels; SRG parameter 
A=1.45 and 1.5 fm _1 . Experimental data (see Ref. [326]. 



Another RGM calculation with the untruncated NCSM approach, but with a SRG transformed 
(A=1.5 fm _1 ) chiral Hamiltonian (I-N3LO) at the two-body level, was carried out for elastic <i- 4 He 
scattering [323]. It is interesting that also a 6 Li ground-state calculation was performed with the RGM 
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approach, where a fragmentation into <i- 4 He and 4 He-(np, T=0) clusters was considered. The virtual 
isoscalar (np)-states were taken into account up to an excitation energy of about 90 MeV, however, the 
RGM result obtained for the 6 Li ground-state energy had not yet converged. In fact it lay 600 keV above 
the result of a six-body NCSM calculation for 6 Li with the same Hamiltonian and where an equivalent 
HO model space was used. As discussed in Ref. [323] the calculation should be further improved by 
taking into account excited states of 4 He and other deuteron excitations not yet considered. Perhaps 
also a fragmentation into 3 H- 3 He could lead to non-negligible effects. For the scattering calculation, i.e. 
for s- and d-wave phase shifts, it was checked that the HO model space (iV max = 12) and the number 
of virtual np states was sufficiently large. However, as already mentioned, unincluded 4 He excitations 
and/or different fragmentation channels could have non- negligible effects. The calculated phase shifts 
for 3 si and 3 <i 3 agree quite well with experimental data, whereas differences were found for 3 di and 3 d 2 . 
In addition, differential cross sections were compared to data at various energies in a range of — 3 - 
12 MeV. A fair agreement of theory and experiment was found. The results obtained are certainly 
very impressive, but it should not be forgotten that in a rigorous ab initio calculation the three-body 
break-up into 4 He-p-n should be included. 

Recently the RGM/NCSM approach was applied to three reactions of astrophysical importance, 
namely to 7 Be(p, 7) 8 B [324] and to the fusion reactions 3 H(rf, n) 4 He (see Fig. 23) and 3 He(ci,p) 4 He [325]. 
The latter two are also of interest for future energy generation on earth. In both cases a SRG evolved 
chiral Hamiltonian (I-N3LO), at the two-body level, was used, hence SRG evolved 3NF as well as bare 
3NF were not considered. In Ref. [324] the SRG parameter A was chosen within what the authors 
called a natural range with values from 1.73 to 1.95 fm _1 . With this choice the 8 B separation energies 
could be well reproduced. On the contrary for the fusion reactions A = 1.5 fm _1 was taken in order to 
better describe the Q values of both reactions. It is evident that such reaction dependent choices of A 
are against the spirit of a true ab initio calculation. This is also admitted by the authors, who state 
that it would be desirable to include a SRG transformation at the three-body level and a consideration 
of a bare 3NF. Then, in addition, it should be either shown that the results do not depend on A or 
alternatively one may define a new potential model choosing a fixed value for A. In spite of all that we 
consider these calculations to represent a considerable progress. 

For the 7 Be(p, 7) 8 -B reaction it was necessary to use the NCSM truncation scheme in order to 
have sufficiently convergent results. The 8 B ground and continuum states were described as (1+7) 
fragmentations, where, obviously, the proton is the single particle, while the 7 Be compound system was 
described by the ground state and four additional excited states. For the radiative capture cross section 
only the dominant El transitions were considered and were calculated by using the unretarded dipole 
approximation (see Section 3.2.4). The resulting S'-factor, 517(0), was found to be approximately 19.4 
eV b, which was consistent with the latest evaluation of 20.8 ± 0.7(expt) ± 1.4(theory) eV b [262]. 

For the calculation of the fusion reactions, first the phase shifts for elastic scattering of n- 4 He, d- 
3 H, and <i- 3 He were determined by including channel coupling, something which was not done in the 
previous NCSM studies of n- 4 He scattering. The effect of the channel coupling on the phase shifts 
remained quite small, but in some cases was not negligible. Also virtual (np, T = 0) states were taken 
into account, which led to sizable effects in the 4 S3/2 phase shifts of the <i- 3 H and d- 3 He channels. The 
results for the S-factors of the fusion reactions showed that the virtual np states led to important effects. 
A rather good agreement with experimental data was obtained in the case of 3 He(d, p) 4 H (approximately 
correct peak position, slight overestimation of peak height), whereas a slightly wrong peak position was 
found for 3 H(<i, n) 4 H. Interesting is the fact, that only a slight change of the SRG parameter A from 
1.5 to 1.45 fm -1 leads to quite a good description of the S'-factor for the 3 H(<i, n) 4 H reaction, but 
presumably would lead to a stronger disagreement in the 3 He(<i, p) 4 H channel. This demonstrates again 
that it would be desirable to perform a more complete SRG transformation and to include a bare 3NF, 
although one has to consider that such calculations are quite demanding. 

Another reaction of astrophysical relevance, 3 He(o;, 7) 7 Be, was calculated by Neff in unretarded 
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Figure 24: Total photoabsorption cross sections of 6 He and 7 Li calculated in Refs. [329] and [330], re- 
spectively, with the LIT/EIHH method and the AV4' [331] potential. Experimental data from Ref. [332]. 



dipole approximation using the FMD approach [327]. As a NN interaction a UCOM version of the AV18 
was taken. Quite a good agreement with experimental data was obtained, both for 3 He - 4 He scattering 
phase shifts and the 3 He(o;,7) 7 Be cross section. The reaction with the mirror nuclei, 3 H(o;,7) 7 Li, was 
calculated as well. The shape of the experimental cross section was reproduced, but the absolute 
normalization was about 15% above the most recent experimental data. However, as already pointed 
out in Section 2.3.5 the FMD technique still needs to be put on safer grounds since benchmark tests 
with other ab initio methods are missing. It would be desirable to have such tests. 

Coming to the end of this section we want to devote some space to ab initio calculations which 
reach further into the many-body continuum. As already pointed out, such calculations can be carried 
out with the LIT method. Unfortunately, LIT calculations have not yet been performed for A > 4 
nuclei with realistic nuclear Hamiltonians although we hope that this will soon be possible at least for 
A = 6, 7. Therefore we want to discuss some examples with simpler forces. In Fig. 24 we show the total 
photoabsorption cross sections of 6 He [328, 329] and 7 Li [330] calculated by Bacca et al. with the central 
NN potential AV4' and using the EIHH method. One notes that the 6 He cross section exhibits a very 
interesting structure. In fact two separate peaks were found. One may interpret the low-energy peak 
as a response due to the relative motion of the two halo neutrons with respect to the a core. Similarly 
one may interpret the peak at higher energies in terms of the classical giant resonance picture, i.e. a 
collective response of protons with respect to the neutrons. For the latter case a break-up of the a core 
is necessary and thus the peak is located at higher energies. A similar structure of the cross section was 
not found for the 7 Li total photoabsorption cross section where only a giant resonance peak is present. 
However, it should be noted that whereas good quality experimental data exist for the 7 Li case the 
same is not true for the 6 He case. The theoretical LIT cross section describes the 7 Li data fairly well 
far up into the many-body continuum, where many different break-up channels should contribute to 
the total cross section. In fact with the LIT method the effect of all different channels is rigorously 
included within an ab initio approach. 



4 Summary and Outlook 



We have presented an overview of modern ab initio methods in few-nucleon systems with A > 4. Our 
review shows that great progress has been made in the past, and in particular in the last two decades. 
As discussed in Section 2, a number of strategies have been followed for solving the nuclear many-body 
problem. Some of the methods have already been introduced rather long ago (see Section 3.1), while 
other approaches have appeared on the scene only rather recently (see Section 3.2). 

As pointed out in the introduction we have qualified a method to be an ab initio method if it allows 
a solution of the nuclear many-body problem without any uncontrolled approximation. Controlled 
expansions, however, can be increasingly improved such that eventually a converged result is obtained. 
Such results are especially relevant for extracting reliable information about one of the fundamentals of 
nuclear physics, i.e. the nuclear force. 

Precise ab initio results from different ab initio methods have to agree within a reasonable error 
range and hence they become benchmark results. In our report we have endeavored to show them 
whenever possible. Benchmark results are important and often represent milestones in the development 
of a field. 

We have divided the discussion of results in four parts. In the first part (Section 3.1) we have 
described calculations for the A = A system from an earlier period (1970- 1990). It is remarkable that 
already in 1982 a first milestone was set with benchmark results for the 4 He binding energy with a 
semi-realistic NN potential (see Table 8). The second part of our discussion (Section 3.2) illustrates 
the more modern period of A = 4 calculations. We place their initiation at the end of the 1980s when 
a GFMC calculation with precise 4 He ab initio results for a realistic NN potential has been published. 
Somewhat later a similar success was achieved by the FY approach. In the following years a great 
innovation phase had its beginnings. Various groups had begun to tackle the 4 He ground-state problem 
with different methods. This culminated in another milestone ten years later - a benchmark calculation 
for the a particle with a realistic NN potential, where seven different groups using seven different 
methods participated (see Tables 10 and 11). 

In Section 3.2 we reported on calculations for the 4N continuum, discussing results for hadronic and 
perturbation induced reactions. These calculations are more complicated than for the 4N ground state, 
thus accounting for the much later appearance of the benchmarks. In fact, only very recently have three 
groups succeeded in doing so, by considering elastic (3+1) scattering below the break-up threshold (see 
Section 3.2.3). The comparison to experiment, however, is still problematic in some cases, even when 
modern realistic nuclear forces are used (see for example the total n- 3 H elastic scattering cross section 
in Fig. 10 or the proton A y in elastic p- 3 He scattering in Fig. 12). These discrepancies call for an 
improvement of the Hamiltonian and need further investigation. 

For higher energies, where not only two- but also more-body channels are open, calculations are 
more complicated and exact computations of such continuum states do not yet exist (with exception 
of the very recent Ref. [153]). However, ab initio calculations of various perturbation-induced reaction 
cross sections, among them the very much discussed 4 He photoabsorption (see Fig. 16 and corresponding 
discussion), have been carried out in this energy range with the LIT method. We think that the power 
of this method, which bypasses the problem of calculating continuum states and aims at computing 
observables directly using bound-state techniques, has not yet been fully explored. In particular the 
application to the LIT of bound-state methods other than the HH or NCSM could be investigated. 
Also the complex scaling approach, which has some similarities to the LIT method, should be further 
examined. 

The energy spectra and other observables of A > 4 systems have been discussed in sections 3.3 and 
3.4, respectively. The calculation of the spectra, pioneered by GFMC has led to a more detailed and 
better description of the 3NF. Improved phenomenological 3NF have allowed a rather high precision 
description of the spectra of 6 < A < 10 nuclei (see also Table 15). We have also discussed quite 
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a number of NCSM studies. Precise ab initio results have been mainly obtained for A < 6 nuclei, 
while for A > 6 nuclei, particularly in the presence of a 3NF, the model space has often not been 
sufficiently large to reach converged results. However, in recent years the NCSM convergence has been 
significantly improved by a controlled truncation of basis states. Also the introduction of SRG evolved 
Hamiltonians has allowed considerable acceleration towards convergence, although still not enough to 
reach full convergence for the ground states of 12 C and 16 in presence of a 3NF. Without a 3NF the 
NCSM ground-state energy of 16 agrees fairly well with CC and UMOA results. Evidently, the CC 
and UMOA techniques are very powerful as well, even if up to now they have been restricted mainly to 
closed shell nuclei. On the other hand, very recently successful CC attempts have been made to work 
in open-shell nuclei. For further observables of A > 4 nuclei (radii, magnetic moments, electromagnetic 
transition strengths, scattering observables, S'-factors, total photoabsorption cross sections) there are a 
wealth of calculations. In 3.4 we have limited ourself to a few selected examples of which we want to 
mention here only the great progress of scattering calculations with the RGM/NCSM approach. 

Considering the progress that has been made in few-nucleon physics over the years, and in particular 
in the last two decades, we think that the future prospects are very promising and that further important 
progress will be made in the years to come. In fact, there still remains a great deal of work to be 
done. For example, it would be desirable to have benchmark tests with realistic Hamiltonians for 
A = 4 continuum calculations beyond the three-body break-up threshold, as well as for ground- and 
continuum-state quantities of A > 4 nuclear systems. Another issue to be addressed concerns SRG 
transformations. In the last decade it has become evident that they play a relevant role in present-day 
few-nucleon physics. One important aspect in such transformations is that they are applied only at the 
two- or three-body level, but the resulting Hamiltonian is used in the A-bodj calculation. Higher order 
terms may be negligible and then the transformation is unitary on the A-body level (independence of 
results on flow parameter), otherwise one has induced many-body forces which do not correspond to 
those of the unitary transformed A-body Hamiltonian. However, a three-body level SRG Hamiltonian 
of a bare two- plus three-nucleon force can be as legitimate as the original one, provided that the flow 
parameter is adjusted once for all and then kept fixed. Otherwise the predictive power of the theory 
would be lost. 

For the advancement of the field important contributions must also come from experiment. While 
there already exists a wealth of data on the bound-state properties of stable nuclei what is needed are 
precise measurements of observables in low-energy reactions. Particularly needed are those observables 
arising from electromagnetic probes (e.g. electron scattering observables at low energy and low mo- 
mentum transfer) since, as pointed out in Section 3.2.2, they offer alternative insights into the nuclear 
dynamics. In addition, experiments with beams of unstable nuclei such as 6 He or other light neutron 
rich nuclei, will be an important source of complementary information. Present and future facilities 
with radioactive ion beams like GANIL /SPIRAL2, GSI/FAIR and RIKEN/RIBF are expected to play 
a key role in the determination of such observables. 

Concluding this review we express our confidence that the field of ab initio calculations will continue 
to grow in importance in future nuclear physics. 
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